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 1861 by gezelter, Tue Apr 9 19:45:54 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 1435 | Line 1445 | namespace OpenMD {
1445      RealType Mc = 0.0;
1446      Mat3x3d Ic(0.0);
1447      RealType Kc = 0.0;
1448 +
1449 +    // Constraints can be on only the linear or angular momentum, but
1450 +    // not both.  Usually, the user will specify which they want, but
1451 +    // in case they don't, the use of periodic boundaries should make
1452 +    // the choice for us.
1453 +    bool doLinearPart = false;
1454 +    bool doAngularPart = false;
1455 +
1456 +    switch (rnemdFluxType_) {
1457 +    case rnemdPx:
1458 +    case rnemdPy:
1459 +    case rnemdPz:
1460 +    case rnemdPvector:
1461 +    case rnemdKePx:
1462 +    case rnemdKePy:
1463 +    case rnemdKePvector:
1464 +      doLinearPart = true;
1465 +      break;
1466 +    case rnemdLx:
1467 +    case rnemdLy:
1468 +    case rnemdLz:
1469 +    case rnemdLvector:
1470 +    case rnemdKeLx:
1471 +    case rnemdKeLy:
1472 +    case rnemdKeLz:
1473 +    case rnemdKeLvector:
1474 +      doAngularPart = true;
1475 +      break;
1476 +    case rnemdKE:
1477 +    case rnemdRotKE:
1478 +    case rnemdFullKE:
1479 +    default:
1480 +      if (usePeriodicBoundaryConditions_)
1481 +        doLinearPart = true;
1482 +      else
1483 +        doAngularPart = true;
1484 +      break;
1485 +    }
1486      
1487      for (sd = smanA.beginSelected(selei); sd != NULL;
1488           sd = smanA.nextSelected(selei)) {
# Line 1543 | Line 1591 | namespace OpenMD {
1591                                MPI::REALTYPE, MPI::SUM);
1592   #endif
1593      
1594 +
1595 +    Vector3d ac, acrec, bc, bcrec;
1596 +    Vector3d ah, ahrec, bh, bhrec;
1597 +    RealType cNumerator, cDenominator;
1598 +    RealType hNumerator, hDenominator;
1599 +
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 >      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 >        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 >            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 >              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 1596 | Line 1669 | namespace OpenMD {
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 1609 | Line 1685 | namespace OpenMD {
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 1751 | Line 1830 | namespace OpenMD {
1830    void RNEMD::collectData() {
1831      if (!doRNEMD_) return;
1832      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1833 <
1833 >    
1834      // collectData can be called more frequently than the doRNEMD, so use the
1835      // computed area from the last exchange time:
1836 <
1837 <    areaAccumulator_->add(getDividingArea());
1836 >    RealType area = getDividingArea();
1837 >    areaAccumulator_->add(area);
1838      Mat3x3d hmat = currentSnap_->getHmat();
1839      seleMan_.setSelectionSet(evaluator_.evaluate());
1840  
# Line 1814 | Line 1893 | namespace OpenMD {
1893        Vector3d rPos = sd->getPos() - coordinateOrigin_;
1894        Vector3d aVel = cross(rPos, vel);
1895        
1896 <      if (binNo < nBins_)  {
1896 >      if (binNo >= 0 && binNo < nBins_)  {
1897          binCount[binNo]++;
1898          binMass[binNo] += mass;
1899          binPx[binNo] += mass*vel.x();
# Line 1890 | Line 1969 | namespace OpenMD {
1969        vel.x() = binPx[i] / binMass[i];
1970        vel.y() = binPy[i] / binMass[i];
1971        vel.z() = binPz[i] / binMass[i];
1972 <      aVel.x() = binOmegax[i];
1973 <      aVel.y() = binOmegay[i];
1974 <      aVel.z() = binOmegaz[i];
1972 >      aVel.x() = binOmegax[i] / binCount[i];
1973 >      aVel.y() = binOmegay[i] / binCount[i];
1974 >      aVel.z() = binOmegaz[i] / binCount[i];
1975  
1976        if (binCount[i] > 0) {
1977          // only add values if there are things to add
# Line 1914 | Line 1993 | namespace OpenMD {
1993              case VELOCITY:
1994                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1995                break;
1996 <            case ANGULARVELOCITY:
1996 >            case ANGULARVELOCITY:  
1997                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
1998                break;
1999              case DENSITY:
# Line 2069 | Line 2148 | namespace OpenMD {
2148            if (outputMask_[i]) {
2149              if (data_[i].dataType == "RealType")
2150                writeReal(i,j);
2151 <            else if (data_[i].dataType == "Vector3d")
2151 >            else if (data_[i].dataType == "Vector3d")
2152                writeVector(i,j);
2153              else {
2154                sprintf( painCave.errMsg,
# Line 2126 | Line 2205 | namespace OpenMD {
2205      RealType s;
2206      int count;
2207      
2208 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2208 >    count = data_[index].accumulator[bin]->count();
2209      if (count == 0) return;
2210      
2211      dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s);
# Line 2149 | Line 2228 | namespace OpenMD {
2228      Vector3d s;
2229      int count;
2230      
2231 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2231 >    count = data_[index].accumulator[bin]->count();
2232 >
2233      if (count == 0) return;
2234  
2235      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
# Line 2173 | Line 2253 | namespace OpenMD {
2253      RealType s;
2254      int count;
2255      
2256 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2256 >    count = data_[index].accumulator[bin]->count();
2257      if (count == 0) return;
2258      
2259      dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
# Line 2196 | Line 2276 | namespace OpenMD {
2276      Vector3d s;
2277      int count;
2278      
2279 <    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count();
2279 >    count = data_[index].accumulator[bin]->count();
2280      if (count == 0) return;
2281  
2282      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);

Diff Legend

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