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 1867 by gezelter, Mon Apr 29 17:53:48 2013 UTC

# Line 419 | Line 419 | namespace OpenMD {
419          velocity.accumulator.push_back( new VectorAccumulator() );
420        data_[VELOCITY] = velocity;
421        outputMap_["VELOCITY"] = VELOCITY;
422 +
423 +      OutputData angularVelocity;
424 +      angularVelocity.units = "angstroms^2/fs";
425 +      angularVelocity.title =  "AngularVelocity";  
426 +      angularVelocity.dataType = "Vector3d";
427 +      angularVelocity.accumulator.reserve(nBins_);
428 +      for (int i = 0; i < nBins_; i++)
429 +        angularVelocity.accumulator.push_back( new VectorAccumulator() );
430 +      data_[ANGULARVELOCITY] = angularVelocity;
431 +      outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
432  
433        OutputData density;
434        density.units =  "g cm^-3";
# Line 546 | Line 556 | namespace OpenMD {
556              slabWidth_ = hmat(2,2) / 10.0;
557          
558            if (hasSlabBCenter)
559 <            slabBCenter_ = rnemdParams->getSlabACenter();
559 >            slabBCenter_ = rnemdParams->getSlabBCenter();
560            else
561              slabBCenter_ = hmat(2,2) / 2.0;
562          
# Line 567 | Line 577 | namespace OpenMD {
577          }
578        }
579      }
580 +
581      // object evaluator:
582      evaluator_.loadScriptString(rnemdObjectSelection_);
583      seleMan_.setSelectionSet(evaluator_.evaluate());
573 –    
584      evaluatorA_.loadScriptString(selectionA_);
585      evaluatorB_.loadScriptString(selectionB_);
576 –    
586      seleManA_.setSelectionSet(evaluatorA_.evaluate());
587      seleManB_.setSelectionSet(evaluatorB_.evaluate());
579 –    
588      commonA_ = seleManA_ & seleMan_;
589 <    commonB_ = seleManB_ & seleMan_;    
589 >    commonB_ = seleManB_ & seleMan_;  
590    }
591    
592      
# Line 1435 | Line 1443 | namespace OpenMD {
1443      RealType Mc = 0.0;
1444      Mat3x3d Ic(0.0);
1445      RealType Kc = 0.0;
1446 +
1447 +    // Constraints can be on only the linear or angular momentum, but
1448 +    // not both.  Usually, the user will specify which they want, but
1449 +    // in case they don't, the use of periodic boundaries should make
1450 +    // the choice for us.
1451 +    bool doLinearPart = false;
1452 +    bool doAngularPart = false;
1453 +
1454 +    switch (rnemdFluxType_) {
1455 +    case rnemdPx:
1456 +    case rnemdPy:
1457 +    case rnemdPz:
1458 +    case rnemdPvector:
1459 +    case rnemdKePx:
1460 +    case rnemdKePy:
1461 +    case rnemdKePvector:
1462 +      doLinearPart = true;
1463 +      break;
1464 +    case rnemdLx:
1465 +    case rnemdLy:
1466 +    case rnemdLz:
1467 +    case rnemdLvector:
1468 +    case rnemdKeLx:
1469 +    case rnemdKeLy:
1470 +    case rnemdKeLz:
1471 +    case rnemdKeLvector:
1472 +      doAngularPart = true;
1473 +      break;
1474 +    case rnemdKE:
1475 +    case rnemdRotKE:
1476 +    case rnemdFullKE:
1477 +    default:
1478 +      if (usePeriodicBoundaryConditions_)
1479 +        doLinearPart = true;
1480 +      else
1481 +        doAngularPart = true;
1482 +      break;
1483 +    }
1484      
1485      for (sd = smanA.beginSelected(selei); sd != NULL;
1486           sd = smanA.nextSelected(selei)) {
# Line 1543 | Line 1589 | namespace OpenMD {
1589                                MPI::REALTYPE, MPI::SUM);
1590   #endif
1591      
1592 +
1593 +    Vector3d ac, acrec, bc, bcrec;
1594 +    Vector3d ah, ahrec, bh, bhrec;
1595 +    RealType cNumerator, cDenominator;
1596 +    RealType hNumerator, hDenominator;
1597 +
1598 +
1599      bool successfulExchange = false;
1600      if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty
1601        Vector3d vc = Pc / Mc;
1602 <      Vector3d ac = -momentumTarget_ / Mc + vc;
1603 <      Vector3d acrec = -momentumTarget_ / Mc;
1602 >      ac = -momentumTarget_ / Mc + vc;
1603 >      acrec = -momentumTarget_ / Mc;
1604        
1605        // We now need the inverse of the inertia tensor to calculate the
1606        // angular velocity of the cold slab;
1607        Mat3x3d Ici = Ic.inverse();
1608        Vector3d omegac = Ici * Lc;
1609 <      Vector3d bc  = -(Ici * angularMomentumTarget_) + omegac;
1610 <      Vector3d bcrec = bc - omegac;
1609 >      bc  = -(Ici * angularMomentumTarget_) + omegac;
1610 >      bcrec = bc - omegac;
1611        
1612 <      RealType cNumerator = Kc - kineticTarget_
1613 <        - 0.5 * Mc * ac.lengthSquare() - 0.5 * ( dot(bc, Ic * bc));
1612 >      cNumerator = Kc - kineticTarget_;
1613 >      if (doLinearPart)
1614 >        cNumerator -= 0.5 * Mc * ac.lengthSquare();
1615 >      
1616 >      if (doAngularPart)
1617 >        cNumerator -= 0.5 * ( dot(bc, Ic * bc));
1618 >
1619        if (cNumerator > 0.0) {
1620          
1621 <        RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare()
1622 <          - 0.5*(dot(omegac, Ic * omegac));
1621 >        cDenominator = Kc;
1622 >
1623 >        if (doLinearPart)
1624 >          cDenominator -= 0.5 * Mc * vc.lengthSquare();
1625 >
1626 >        if (doAngularPart)
1627 >          cDenominator -= 0.5*(dot(omegac, Ic * omegac));
1628          
1629          if (cDenominator > 0.0) {
1630            RealType c = sqrt(cNumerator / cDenominator);
1631            if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients
1632              
1633              Vector3d vh = Ph / Mh;
1634 <            Vector3d ah = momentumTarget_ / Mh + vh;
1635 <            Vector3d ahrec = momentumTarget_ / Mh;
1634 >            ah = momentumTarget_ / Mh + vh;
1635 >            ahrec = momentumTarget_ / Mh;
1636              
1637              // We now need the inverse of the inertia tensor to
1638              // calculate the angular velocity of the hot slab;
1639              Mat3x3d Ihi = Ih.inverse();
1640              Vector3d omegah = Ihi * Lh;
1641 <            Vector3d bh  = (Ihi * angularMomentumTarget_) + omegah;
1642 <            Vector3d bhrec = bh - omegah;
1641 >            bh  = (Ihi * angularMomentumTarget_) + omegah;
1642 >            bhrec = bh - omegah;
1643              
1644 <            RealType hNumerator = Kh + kineticTarget_
1645 <              - 0.5 * Mh * ah.lengthSquare() - 0.5 * ( dot(bh, Ih * bh));;
1644 >            hNumerator = Kh + kineticTarget_;
1645 >            if (doLinearPart)
1646 >              hNumerator -= 0.5 * Mh * ah.lengthSquare();
1647 >            
1648 >            if (doAngularPart)
1649 >              hNumerator -= 0.5 * ( dot(bh, Ih * bh));
1650 >              
1651              if (hNumerator > 0.0) {
1652                
1653 <              RealType hDenominator = Kh - 0.5 * Mh * vh.lengthSquare()
1654 <                - 0.5*(dot(omegah, Ih * omegah));
1653 >              hDenominator = Kh;
1654 >              if (doLinearPart)
1655 >                hDenominator -= 0.5 * Mh * vh.lengthSquare();
1656 >              if (doAngularPart)
1657 >                hDenominator -= 0.5*(dot(omegah, Ih * omegah));
1658                
1659                if (hDenominator > 0.0) {
1660                  RealType h = sqrt(hNumerator / hDenominator);
# Line 1596 | Line 1667 | namespace OpenMD {
1667                    for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
1668                      //vel = (*sdi)->getVel();
1669                      rPos = (*sdi)->getPos() - coordinateOrigin_;
1670 <                    vel = ((*sdi)->getVel() - vc - cross(omegac, rPos)) * c
1671 <                      + ac + cross(bc, rPos);
1670 >                    if (doLinearPart)
1671 >                      vel = ((*sdi)->getVel() - vc) * c + ac;
1672 >                    if (doAngularPart)
1673 >                      vel = ((*sdi)->getVel() - cross(omegac, rPos)) * c + cross(bc, rPos);
1674 >
1675                      (*sdi)->setVel(vel);
1676                      if (rnemdFluxType_ == rnemdFullKE) {
1677                        if ((*sdi)->isDirectional()) {
# Line 1609 | Line 1683 | namespace OpenMD {
1683                    for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
1684                      //vel = (*sdi)->getVel();
1685                      rPos = (*sdi)->getPos() - coordinateOrigin_;
1686 <                    vel = ((*sdi)->getVel() - vh - cross(omegah, rPos)) * h
1687 <                      + ah + cross(bh, rPos);                  
1686 >                    if (doLinearPart)
1687 >                      vel = ((*sdi)->getVel() - vh) * h + ah;    
1688 >                    if (doAngularPart)
1689 >                      vel = ((*sdi)->getVel() - cross(omegah, rPos)) * h + cross(bh, rPos);    
1690 >
1691                      (*sdi)->setVel(vel);
1692                      if (rnemdFluxType_ == rnemdFullKE) {
1693                        if ((*sdi)->isDirectional()) {
# Line 1653 | Line 1730 | namespace OpenMD {
1730        int isd;
1731        StuntDouble* sd;
1732        vector<StuntDouble*> aSites;
1656 –      ConvexHull* surfaceMeshA = new ConvexHull();
1733        seleManA_.setSelectionSet(evaluatorA_.evaluate());
1734        for (sd = seleManA_.beginSelected(isd); sd != NULL;
1735             sd = seleManA_.nextSelected(isd)) {
1736          aSites.push_back(sd);
1737        }
1738 +      ConvexHull* surfaceMeshA = new ConvexHull();
1739        surfaceMeshA->computeHull(aSites);
1740        areaA = surfaceMeshA->getArea();
1741 +      delete surfaceMeshA;
1742 +
1743      } else {
1744        if (usePeriodicBoundaryConditions_) {
1745          // in periodic boundaries, the surface area is twice the x-y
# Line 1674 | Line 1753 | namespace OpenMD {
1753        }
1754      }
1755  
1756 +
1757 +
1758      if (hasSelectionB_) {
1759        int isd;
1760        StuntDouble* sd;
1761        vector<StuntDouble*> bSites;
1681 –      ConvexHull* surfaceMeshB = new ConvexHull();
1762        seleManB_.setSelectionSet(evaluatorB_.evaluate());
1763        for (sd = seleManB_.beginSelected(isd); sd != NULL;
1764             sd = seleManB_.nextSelected(isd)) {
1765          bSites.push_back(sd);
1766        }
1767 +      ConvexHull* surfaceMeshB = new ConvexHull();    
1768        surfaceMeshB->computeHull(bSites);
1769        areaB = surfaceMeshB->getArea();
1770 +      delete surfaceMeshB;
1771 +
1772      } else {
1773        if (usePeriodicBoundaryConditions_) {
1774          // in periodic boundaries, the surface area is twice the x-y
# Line 1751 | Line 1834 | namespace OpenMD {
1834    void RNEMD::collectData() {
1835      if (!doRNEMD_) return;
1836      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1837 <
1837 >    
1838      // collectData can be called more frequently than the doRNEMD, so use the
1839      // computed area from the last exchange time:
1840 <
1841 <    areaAccumulator_->add(getDividingArea());
1840 >    RealType area = getDividingArea();
1841 >    areaAccumulator_->add(area);
1842      Mat3x3d hmat = currentSnap_->getHmat();
1843      seleMan_.setSelectionSet(evaluator_.evaluate());
1844  
# Line 1814 | Line 1897 | namespace OpenMD {
1897        Vector3d rPos = sd->getPos() - coordinateOrigin_;
1898        Vector3d aVel = cross(rPos, vel);
1899        
1900 <      if (binNo < nBins_)  {
1900 >      if (binNo >= 0 && binNo < nBins_)  {
1901          binCount[binNo]++;
1902          binMass[binNo] += mass;
1903          binPx[binNo] += mass*vel.x();
# Line 1890 | Line 1973 | namespace OpenMD {
1973        vel.x() = binPx[i] / binMass[i];
1974        vel.y() = binPy[i] / binMass[i];
1975        vel.z() = binPz[i] / binMass[i];
1976 <      aVel.x() = binOmegax[i];
1977 <      aVel.y() = binOmegay[i];
1978 <      aVel.z() = binOmegaz[i];
1976 >      aVel.x() = binOmegax[i] / binCount[i];
1977 >      aVel.y() = binOmegay[i] / binCount[i];
1978 >      aVel.z() = binOmegaz[i] / binCount[i];
1979  
1980        if (binCount[i] > 0) {
1981          // only add values if there are things to add
# Line 1914 | Line 1997 | namespace OpenMD {
1997              case VELOCITY:
1998                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1999                break;
2000 <            case ANGULARVELOCITY:
2000 >            case ANGULARVELOCITY:  
2001                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
2002                break;
2003              case DENSITY:
# Line 2069 | Line 2152 | namespace OpenMD {
2152            if (outputMask_[i]) {
2153              if (data_[i].dataType == "RealType")
2154                writeReal(i,j);
2155 <            else if (data_[i].dataType == "Vector3d")
2155 >            else if (data_[i].dataType == "Vector3d")
2156                writeVector(i,j);
2157              else {
2158                sprintf( painCave.errMsg,
# Line 2126 | Line 2209 | namespace OpenMD {
2209      RealType s;
2210      int count;
2211      
2212 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2212 >    count = data_[index].accumulator[bin]->count();
2213      if (count == 0) return;
2214      
2215      dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s);
# Line 2149 | Line 2232 | namespace OpenMD {
2232      Vector3d s;
2233      int count;
2234      
2235 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2235 >    count = data_[index].accumulator[bin]->count();
2236 >
2237      if (count == 0) return;
2238  
2239      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
# Line 2173 | Line 2257 | namespace OpenMD {
2257      RealType s;
2258      int count;
2259      
2260 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2260 >    count = data_[index].accumulator[bin]->count();
2261      if (count == 0) return;
2262      
2263      dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
# Line 2196 | Line 2280 | namespace OpenMD {
2280      Vector3d s;
2281      int count;
2282      
2283 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2283 >    count = data_[index].accumulator[bin]->count();
2284      if (count == 0) return;
2285  
2286      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);

Diff Legend

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