.. _user-guide-rate: Calculating the rate ==================== Everything a path-sampling run does is in service of one number: the rate constant :math:`k_{AB}` at which the system leaves state :math:`A` for state :math:`B`. This page explains how that number is put together -- the identity the methods rest on, where each factor comes from in a |pyretis| run, and how the routes differ between TIS, RETIS, infinite swapping and the (RE)PPTIS family. .. seealso:: :ref:`user-guide-output` explains the files these numbers are computed from and how to read the report; :ref:`user-guide-analyse` is the reference for the application itself. .. contents:: Contents :local: :depth: 2 .. _user-guide-rate-theory: The identity behind every method -------------------------------- A rare transition is rare because the system almost never gets from :math:`A` to :math:`B`. Waiting for it is hopeless; the trick every method here uses is to *factorise* the event into pieces that are each common enough to measure. Put an ordered set of interfaces :math:`\lambda_0 < \lambda_1 < \dots < \lambda_n` between the two states, with :math:`\lambda_0` on the edge of :math:`A` and :math:`\lambda_n` on the edge of :math:`B`. Then .. math:: k_{AB} = f_A \times P_A(\lambda_n | \lambda_0) where * :math:`f_A` is the **initial flux**: how often, per unit time, the system leaves :math:`A` through :math:`\lambda_0` at all. This is a *frequent* event, so plain molecular dynamics measures it. * :math:`P_A(\lambda_n | \lambda_0)` is the **crossing probability**: given that a trajectory has just crossed :math:`\lambda_0` on its way out, the probability it reaches :math:`\lambda_n` before returning to :math:`A`. This is the rare part. The crossing probability is itself factorised, which is the whole point: .. math:: P_A(\lambda_n | \lambda_0) = \prod_{i=0}^{n-1} P_A(\lambda_{i+1} | \lambda_i) Each factor is the probability of getting from one interface to the next -- a *conditional* probability of order 0.1, which ordinary sampling can measure. Their product can be :math:`10^{-7}` or smaller without any single simulation ever having to wait that long. Each ensemble :math:`[i^+]` samples one factor. That is why a rate of :math:`10^{-7}` comes out of runs of a few thousand cycles, and also why the *errors multiply*: see the worked example in :ref:`user-guide-output-worked-example`, where seven factors with 7--34 % relative error each give a rate with 54 %. .. _user-guide-rate-in-code: Where each factor comes from in a run -------------------------------------- The initial flux ~~~~~~~~~~~~~~~~ The flux is measured from the two innermost ensembles, ``[0^-]`` and ``[0^+]``, which between them cover a full round trip in and out of :math:`A`. |pyretis| computes it as the reciprocal of the mean time for that round trip: .. math:: f_A = \frac{1}{(\langle L_{0^-}\rangle - 2 + \langle L_{0^+}\rangle - 2) \, \Delta t} where :math:`\langle L \rangle` is the mean accepted path length in frames and :math:`\Delta t` the time between frames (the engine time step times ``subcycles``). The two subtracted frames are the crossing points shared with the interface, counted once rather than twice. This is the cheap half of the calculation: the flux is usually well determined long before the crossing probability is, which is why a short run reports a sensible :math:`f_A` next to a useless :math:`P_A`. A pure ``task = "md-flux"`` run measures only this, and by a different route: it counts effective crossings of each interface during plain MD and divides by time. There the time in state :math:`A` is cut into windows of :ref:`skipcross ` steps and the flux computed in each, which is what gives the estimate an error bar. The window changes the error, not the flux. The crossing probability ~~~~~~~~~~~~~~~~~~~~~~~~ Each ensemble's own crossing probability curve is a histogram over the highest order parameter its accepted paths reached. Turning the set of curves into one number is called **matching**: each ensemble's curve is scaled so that it continues its neighbour's where they overlap, and the product of the per-ensemble factors is accumulated as the curves are joined. The result is the ``matched-probability`` figure, and its value at the outermost interface is :math:`P_A(\lambda_n|\lambda_0)`. The rate ~~~~~~~~ .. math:: k_{AB} = f_A \times P_A, \qquad \frac{\sigma_k}{k} = \sqrt{\left(\frac{\sigma_f}{f}\right)^2 + \left(\frac{\sigma_P}{P}\right)^2} The relative errors add in quadrature, so the rate is never more precise than its worse factor -- in practice always the crossing probability. .. note:: The rate carries the engine's time unit: a report from an external engine says ``1/gromacs`` or similar. Convert before comparing. .. _user-guide-rate-per-method: How the methods differ ----------------------- All of them sample paths and all of them end at :math:`k = f_A P_A`. They differ in *which* ensembles they define and *how* the per-ensemble probabilities are combined. .. list-table:: :header-rows: 1 :widths: 14 44 42 * - ``task`` - What it samples - How the rate is formed * - ``md-flux`` - Plain MD in :math:`A`, recording crossings of each interface (``cross.txt``). - Gives :math:`f_A` only -- no crossing probability, so no rate. Use it to get the flux cheaply, or to place interfaces. * - ``tis`` - One path ensemble per interface, each sampled independently. - The per-ensemble probabilities are matched and multiplied exactly as above. With no ``[0^-]`` ensemble there is no flux, so a pure TIS run reports a crossing probability and needs the flux from a separate ``md-flux`` run. * - ``retis`` - The same ensembles plus ``[0^-]``, sampled *together*, with swap moves exchanging paths between neighbours. - Flux and crossing probability come from the one run, so it reports the rate directly. This is the usual choice. * - ``infinite_swapping`` - The same physics as RETIS, but instead of accepting or rejecting a swap the coordinator spreads each path's *occupancy* over the ensembles it is valid in. - Identical sampling target, so the same matched estimator applies. A second estimator (WHAM) is also available -- see below. * - ``pptis`` / ``repptis`` - **Partial** paths: each ensemble samples only between neighbouring interfaces rather than all the way back to :math:`A`, which keeps path lengths short in long or diffusive systems. - The local probabilities are assembled into a Markov state model over the interfaces and the global crossing probability is obtained from it (:py:func:`~pyretis.analysis.repptis_msm.global_pcross_msm`), rather than by a plain product. The same MSM also yields mean first-passage times. .. note:: Whichever ``task`` is written in the input, the *sampling* goes through one code path -- the infinite-swapping coordinator. What the task selects is the ensemble layout, the move set and the estimator. .. _user-guide-rate-estimators: Two estimators over the same data ---------------------------------- For a run that produced per-ensemble output, two independent estimators can be applied to it, chosen with the :ref:`method ` keyword: ``matched`` (the default; ``crossing`` is accepted as its old name) The point-matching product described above -- the classic RETIS estimator, with block-error analysis. ``wham`` Weighted-histogram analysis over the same paths, which reweights all ensembles together instead of matching them pairwise. ``both`` Run each and report both. They are different estimators of the *same* physical rate, so on a converged run they agree within their statistical error. Disagreement beyond that is a signal -- usually that the run is not as converged as the error bars suggest. Running ``method = "both"`` is a cheap cross-check, since both read the same files. .. _user-guide-rate-running: Running the analysis -------------------- .. code-block:: bash pyretis analyse -i input.toml The input is the same file the simulation used (or the ``output.toml`` the run wrote, which carries the resolved settings). The analysis reads the per-ensemble output in place and writes ``report/``. The arguments are listed in :ref:`tableappargument_analyse`; the two worth knowing are ``-skipb`` and ``-skipe``, which choose the cycle window without editing the input file: .. code-block:: bash pyretis analyse -i input.toml -skipb 500 # drop equilibration pyretis analyse -i input.toml -skipb 500 -skipe 100 # and a cut-short tail Everything else is set in the :ref:`[analysis] section ` of the input file. The keywords that change the *result* rather than the presentation are: .. list-table:: :header-rows: 1 :widths: 26 74 * - Keyword - Effect on the rate * - ``method`` - Which estimator(s) to use: ``matched``, ``wham`` or ``both``. * - ``skip_initial_cycles`` - Cycles discarded at the start; ``-skipb`` overrides it. * - ``skip_final_cycles`` - Cycles discarded at the end; ``-skipe`` overrides it. * - ``maxblock``, ``blockskip`` - The block lengths used for the error estimate. They do not change the rate, but they do change its quoted error -- and the plateau of the block-error curve is what makes that error trustworthy. * - ``ngrid``, ``bins`` - Resolution of the crossing-probability grid and of the histograms. Too coarse a grid visibly steps the matched curve. A run analysed twice with different windows gives two legitimate answers over two different sets of cycles; the log states which cycles were used, and the report is stamped with the cycle count. .. _user-guide-rate-sanity: Is the number trustworthy? -------------------------- The rate is the most-processed quantity the code produces, so it is the last thing to converge and the first to mislead. Before quoting one: #. Work through :ref:`user-guide-output-converged` -- especially whether the running average has flattened and the block-error curve has reached its plateau. #. Check the flux and the crossing probability *separately*. They fail differently: the flux is usually fine early, so a suspicious rate is nearly always a crossing-probability problem. #. Look at the matched-probability curve for kinks at the interfaces. A kink means neighbouring ensembles disagree where they overlap, and the product across that seam is not to be trusted. #. If in doubt, run ``method = "both"`` and compare the two estimators. |pyretis| ships a validation suite (``examples/validation/``) that runs this whole chain on a one-dimensional double well whose rate is known analytically from Kramers' theory, and checks that the sampled rate lands on it. It is the reference for "does the machinery still give the right answer" -- see that directory's ``README.rst``.