| 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(); |