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

Comparing branches/development/src/rnemd/RNEMD.cpp (file contents):
Revision 1854 by gezelter, Thu Mar 28 20:54:06 2013 UTC vs.
Revision 1875 by gezelter, Fri May 17 14:41:42 2013 UTC

# Line 71 | Line 71 | namespace OpenMD {
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;
# Line 420 | Line 421 | namespace OpenMD {
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";
436        density.title =  "Density";
# Line 546 | Line 557 | namespace OpenMD {
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          
# Line 567 | Line 578 | namespace OpenMD {
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      
# Line 595 | Line 604 | namespace OpenMD {
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) {
# Line 1138 | Line 1151 | namespace OpenMD {
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);
# Line 1149 | Line 1162 | namespace OpenMD {
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;
# Line 1268 | Line 1281 | namespace OpenMD {
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 ) {
# Line 1300 | Line 1313 | namespace OpenMD {
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_) {
# Line 1367 | Line 1380 | namespace OpenMD {
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;
# Line 1378 | Line 1391 | namespace OpenMD {
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;
# Line 1435 | Line 1448 | namespace OpenMD {
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)) {
# Line 1543 | Line 1594 | namespace OpenMD {
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 <              - 0.5 * Mh * ah.lengthSquare() - 0.5 * ( dot(bh, Ih * bh));;
1646 >            RealType hNumerator = Kh + kineticTarget_;
1647 >            if (doLinearPart)
1648 >              hNumerator -= 0.5 * Mh * ah.lengthSquare();
1649 >            
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);
# Line 1593 | Line 1666 | namespace OpenMD {
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()) {
# Line 1606 | Line 1682 | namespace OpenMD {
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()) {
# Line 1653 | Line 1732 | namespace OpenMD {
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 +      ConvexHull* surfaceMeshA = new ConvexHull();
1741        surfaceMeshA->computeHull(aSites);
1742        areaA = surfaceMeshA->getArea();
1743 +      delete surfaceMeshA;
1744 +
1745      } else {
1746        if (usePeriodicBoundaryConditions_) {
1747          // in periodic boundaries, the surface area is twice the x-y
# Line 1674 | Line 1755 | namespace OpenMD {
1755        }
1756      }
1757  
1758 +
1759 +
1760      if (hasSelectionB_) {
1761        int isd;
1762        StuntDouble* sd;
1763        vector<StuntDouble*> bSites;
1681 –      ConvexHull* surfaceMeshB = new ConvexHull();
1764        seleManB_.setSelectionSet(evaluatorB_.evaluate());
1765        for (sd = seleManB_.beginSelected(isd); sd != NULL;
1766             sd = seleManB_.nextSelected(isd)) {
1767          bSites.push_back(sd);
1768        }
1769 +      ConvexHull* surfaceMeshB = new ConvexHull();    
1770        surfaceMeshB->computeHull(bSites);
1771        areaB = surfaceMeshB->getArea();
1772 +      delete surfaceMeshB;
1773 +
1774      } else {
1775        if (usePeriodicBoundaryConditions_) {
1776          // in periodic boundaries, the surface area is twice the x-y
# Line 1751 | Line 1836 | namespace OpenMD {
1836    void RNEMD::collectData() {
1837      if (!doRNEMD_) return;
1838      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1839 <
1839 >    
1840      // collectData can be called more frequently than the doRNEMD, so use the
1841      // computed area from the last exchange time:
1842 <
1843 <    areaAccumulator_->add(getDividingArea());
1842 >    RealType area = getDividingArea();
1843 >    areaAccumulator_->add(area);
1844      Mat3x3d hmat = currentSnap_->getHmat();
1845      seleMan_.setSelectionSet(evaluator_.evaluate());
1846  
# Line 1814 | Line 1899 | namespace OpenMD {
1899        Vector3d rPos = sd->getPos() - coordinateOrigin_;
1900        Vector3d aVel = cross(rPos, vel);
1901        
1902 <      if (binNo < nBins_)  {
1902 >      if (binNo >= 0 && binNo < nBins_)  {
1903          binCount[binNo]++;
1904          binMass[binNo] += mass;
1905          binPx[binNo] += mass*vel.x();
# Line 1890 | Line 1975 | namespace OpenMD {
1975        vel.x() = binPx[i] / binMass[i];
1976        vel.y() = binPy[i] / binMass[i];
1977        vel.z() = binPz[i] / binMass[i];
1978 <      aVel.x() = binOmegax[i];
1979 <      aVel.y() = binOmegay[i];
1980 <      aVel.z() = binOmegaz[i];
1978 >      aVel.x() = binOmegax[i] / binCount[i];
1979 >      aVel.y() = binOmegay[i] / binCount[i];
1980 >      aVel.z() = binOmegaz[i] / binCount[i];
1981  
1982        if (binCount[i] > 0) {
1983          // only add values if there are things to add
# Line 1914 | Line 1999 | namespace OpenMD {
1999              case VELOCITY:
2000                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
2001                break;
2002 <            case ANGULARVELOCITY:
2002 >            case ANGULARVELOCITY:  
2003                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
2004                break;
2005              case DENSITY:
# Line 1925 | Line 2010 | namespace OpenMD {
2010          }
2011        }
2012      }
2013 +    hasData_ = true;
2014    }
2015  
2016    void RNEMD::getStarted() {
# Line 1957 | Line 2043 | namespace OpenMD {
2043    
2044    void RNEMD::writeOutputFile() {
2045      if (!doRNEMD_) return;
2046 +    if (!hasData_) return;
2047      
2048   #ifdef IS_MPI
2049      // If we're the root node, should we print out the results
# Line 1978 | Line 2065 | namespace OpenMD {
2065        RealType time = currentSnap_->getTime();
2066        RealType avgArea;
2067        areaAccumulator_->getAverage(avgArea);
1981 –      RealType Jz = kineticExchange_ / (time * avgArea)
1982 –        / PhysicalConstants::energyConvert;
1983 –      Vector3d JzP = momentumExchange_ / (time * avgArea);      
1984 –      Vector3d JzL = angularMomentumExchange_ / (time * avgArea);      
2068  
2069 +      RealType Jz(0.0);
2070 +      Vector3d JzP(V3Zero);
2071 +      Vector3d JzL(V3Zero);
2072 +      if (time >= info_->getSimParams()->getDt()) {
2073 +        Jz = kineticExchange_ / (time * avgArea)
2074 +          / PhysicalConstants::energyConvert;
2075 +        JzP = momentumExchange_ / (time * avgArea);
2076 +        JzL = angularMomentumExchange_ / (time * avgArea);
2077 +      }
2078 +
2079        rnemdFile_ << "#######################################################\n";
2080        rnemdFile_ << "# RNEMD {\n";
2081  
# Line 2069 | Line 2162 | namespace OpenMD {
2162            if (outputMask_[i]) {
2163              if (data_[i].dataType == "RealType")
2164                writeReal(i,j);
2165 <            else if (data_[i].dataType == "Vector3d")
2165 >            else if (data_[i].dataType == "Vector3d")
2166                writeVector(i,j);
2167              else {
2168                sprintf( painCave.errMsg,
# Line 2126 | Line 2219 | namespace OpenMD {
2219      RealType s;
2220      int count;
2221      
2222 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2222 >    count = data_[index].accumulator[bin]->count();
2223      if (count == 0) return;
2224      
2225      dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s);
# Line 2135 | Line 2228 | namespace OpenMD {
2228        rnemdFile_ << "\t" << s;
2229      } else{
2230        sprintf( painCave.errMsg,
2231 <               "RNEMD detected a numerical error writing: %s for bin %d",
2231 >               "RNEMD detected a numerical error writing: %s for bin %u",
2232                 data_[index].title.c_str(), bin);
2233        painCave.isFatal = 1;
2234        simError();
# Line 2149 | Line 2242 | namespace OpenMD {
2242      Vector3d s;
2243      int count;
2244      
2245 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2245 >    count = data_[index].accumulator[bin]->count();
2246 >
2247      if (count == 0) return;
2248  
2249      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
# Line 2157 | Line 2251 | namespace OpenMD {
2251          isinf(s[1]) || isnan(s[1]) ||
2252          isinf(s[2]) || isnan(s[2]) ) {      
2253        sprintf( painCave.errMsg,
2254 <               "RNEMD detected a numerical error writing: %s for bin %d",
2254 >               "RNEMD detected a numerical error writing: %s for bin %u",
2255                 data_[index].title.c_str(), bin);
2256        painCave.isFatal = 1;
2257        simError();
# Line 2173 | Line 2267 | namespace OpenMD {
2267      RealType s;
2268      int count;
2269      
2270 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2270 >    count = data_[index].accumulator[bin]->count();
2271      if (count == 0) return;
2272      
2273      dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
# Line 2182 | Line 2276 | namespace OpenMD {
2276        rnemdFile_ << "\t" << s;
2277      } else{
2278        sprintf( painCave.errMsg,
2279 <               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
2279 >               "RNEMD detected a numerical error writing: %s std. dev. for bin %u",
2280                 data_[index].title.c_str(), bin);
2281        painCave.isFatal = 1;
2282        simError();
# Line 2196 | Line 2290 | namespace OpenMD {
2290      Vector3d s;
2291      int count;
2292      
2293 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2293 >    count = data_[index].accumulator[bin]->count();
2294      if (count == 0) return;
2295  
2296      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
# Line 2204 | Line 2298 | namespace OpenMD {
2298          isinf(s[1]) || isnan(s[1]) ||
2299          isinf(s[2]) || isnan(s[2]) ) {      
2300        sprintf( painCave.errMsg,
2301 <               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
2301 >               "RNEMD detected a numerical error writing: %s std. dev. for bin %u",
2302                 data_[index].title.c_str(), bin);
2303        painCave.isFatal = 1;
2304        simError();

Diff Legend

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