Adaptive biasing with AWH
-------------------------
-The accelerated weight histogram method
+The accelerated weight histogram method :ref:`185 <refLidmar2012>`
:ref:`137 <reflindahl2014accelerated>` calculates the PMF along a reaction coordinate by adding
an adaptively determined biasing potential. AWH flattens free energy
barriers along the reaction coordinate by applying a history-dependent
^^^^^^^^^^^^^^^^^^^^
Rather than biasing the reaction coordinate :math:`\xi(x)` directly, AWH
-acts on a *reference coordinate* :math:`\lambda`. The reaction
-coordinate :math:`\xi(x)` is coupled to :math:`\lambda` with a harmonic
-potential
-
-.. math:: Q(\xi,\lambda) = \frac{1}{2} \beta k (\xi - \lambda)^2,
- :label: eqnawhbasic
-
-so that for large force constants :math:`k`,
-:math:`\xi \approx \lambda`. Note the use of dimensionless energies for
-compatibility with previously published work. Units of energy are
-obtained by multiplication with :math:`k_BT=1/\beta`. In the simulation,
-:math:`\lambda` samples the user-defined sampling interval :math:`I`.
-For a multidimensional reaction coordinate :math:`\xi`, the sampling
-interval is the Cartesian product :math:`I=\Pi_{\mu} I_{\mu}` (a rectangular
-domain). The connection between atom coordinates and :math:`\lambda` is
-established through the extended ensemble \ :ref:`68 <refLyubartsev1992>`,
+acts on a *reference coordinate* :math:`\lambda`. The fundamentals of the
+method is based on the connection between atom coordinates and :math:`\lambda` and
+is established through the extended ensemble \ :ref:`68 <refLyubartsev1992>`,
.. math:: P(x,\lambda) = \frac{1}{\mathcal{Z}}e^{g(\lambda) - Q(\xi(x),\lambda) - V(x)},
:label: eqawhpxlambda
.. math:: F(\lambda) = -\ln \int e^{- Q(\xi(x),\lambda) - V(x)} dx.
:label: eqawhflambda
+The reaction coordinate :math:`\xi(x)` is commonly coupled to
+:math:`\lambda` with a harmonic potential
+
+.. math:: Q(\xi,\lambda) = \frac{1}{2} \beta k (\xi - \lambda)^2,
+ :label: eqnawhbasic
+
+so that for large force constants :math:`k`,
+:math:`\xi \approx \lambda`. Note the use of dimensionless energies for
+compatibility with previously published work. Units of energy are
+obtained by multiplication with :math:`k_BT=1/\beta`. In the simulation,
+:math:`\lambda` samples the user-defined sampling interval :math:`I`.
+
Being the convolution of the PMF with the Gaussian defined by the
harmonic potential, :math:`F(\lambda)` is a smoothened version of the
PMF. :eq:`Eq. %s <eqawhplambda>` shows that in order to obtain
:math:`F(\lambda)` and simultaneously converges :math:`P(\lambda)`
toward :math:`\rho(\lambda)`.
+It is also possible to directly control the :math:`\lambda` state
+of, e.g., alchemical free energy perturbations :ref:`187 <reflundborg2021>`. In that case there is no harmonic
+potential and :math:`\lambda` changes in discrete steps along the reaction coordinate
+depending on the biased free energy difference between the :math:`\lambda` states.
+N.b., it is not yet possible to use AWH in combination with perturbed masses or
+constraints.
+
+For a multidimensional reaction coordinate :math:`\xi`, the sampling
+interval is the Cartesian product :math:`I=\Pi_{\mu} I_{\mu}` (a rectangular
+domain).
+
The free energy update
^^^^^^^^^^^^^^^^^^^^^^
Secondly, the normalization of the histogram
:math:`N_n=\sum_\lambda W_n(\lambda)`, determines the update size
:math:`| \Delta F(\lambda) |`. For instance, for a single sample
-:math:`\omega(\lambda|x)`, the shape of the update is approximately a
+:math:`\omega(\lambda|x)`, and using a harmonic potential
+(:see :eq:`Eq. %s <eqnawhbasic>`),
+the shape of the update is approximately a
Gaussian function of width :math:`\sigma=1/\sqrt{\beta k}` and height
:math:`\propto 1/N_n` :ref:`137 <reflindahl2014accelerated>`,
.. math:: | \Delta F_n(\lambda) | \propto \frac{1}{N_n} e^{-\frac{1}{2} \beta k (\xi(x) - \lambda)^2}.
:label: eqawhdfsize
-Therefore, as samples accumulate in :math:`W(\lambda)` and :math:`N_n`
-grows, the updates get smaller, allowing for the free energy to
+When directly controlling the lambda state of the system, the shape of
+the update is instead
+
+.. math:: | \Delta F_n(\lambda) | \propto \frac{1}{N_n} P_n(\lambda | x).
+ :label: eqawhdfsizelambda
+
+Therefore, in both cases, as samples accumulate in :math:`W(\lambda)` and
+:math:`N_n` grows, the updates get smaller, allowing for the free energy to
converge.
Note that quantity of interest to the user is not :math:`F(\lambda)` but
the AWH algorithm, only what force and potential AWH returns to the MD
engine.
+Along a bias dimension directly controlling the
+:math:`\lambda` state of the system, such as when controlling free energy
+perturbations, the Monte-Carlo sampling alternative is always used, even if
+a convolved bias potential is chosen to be used along the other dimensions
+(if there are more than one).
+
+
.. _fig-awhbiasevolution1:
.. figure:: plots/awh-traj.*
a too rapid decay of the free energy update size. This motivates
splitting the simulation into an *initial stage* where the weight
histogram grows according to a more restrictive and robust protocol, and
-a *final stage* where the the weight histogram grows linearly at the
+a *final stage* where the weight histogram grows linearly at the
sampling rate (:eq:`Eq. %s <eqawhwupdate>`). The AWH initial
stage takes inspiration from the well-known Wang-Landau algorithm \ :ref:`138 <refwang2001efficient>`,
although there are differences in the details.
considered covered when at least one walker has independently traversed
the sampling interval.
+In practice biases are shared by setting :mdp:`awh-share-multisim` to true
+and :mdp:`awh1-share-group` (for bias 1) to a non-zero value. Here, bias 1
+will be shared between simulations that have the same share group value.
+Sharing can be different for bias 1, 2, etc. (although there are
+few use cases where this is useful). Technically there are no restrictions
+on sharing, apart from that biases that are shared need to have the same
+number of grid points and the update intervals should match.
+Note that biases can not be shared within a simulation.
+The latter could be useful, especially for multimeric proteins, but this
+is more difficult to implement.
+
.. _awhreweight:
Reweighting and combining biased data
is the deviation of the force. The factors :math:`\omega(\lambda|x(t))`,
see :eq:`Eq %s <eqawhomega>`, reweight the samples.
:math:`\eta_{\mu\nu}(\lambda)` is a friction
-tensor \ :ref:`144 <refsivak2012thermodynamic>`. Its matrix elements are inversely proportional to local
+tensor :ref:`186 <reflindahl2018>` and :ref:`144 <refsivak2012thermodynamic>`. Its matrix elements are inversely proportional to local
diffusion coefficients. A measure of sampling (in)efficiency at each
:math:`\lambda` is given by
estimate used by AWH to initialize :math:`N` in terms of more meaningful
quantities
-.. math:: \frac{1}{N_0} = \frac{1}{N_0(\varepsilon_0, D)} \sim D\varepsilon_0^2.
+.. math:: \frac{1}{N_0} = \frac{1}{N_0(\varepsilon_0, D)} = \frac{1}{\Delta
+ t_\mathrm{sample}} \max_d \frac{L_d^2}{2D_d} \varepsilon_0^2
:label: eqawhn0
-Essentially, this tells us that a slower system (small :math:`D`)
+where :math:`L_d` is the length of the interval and :math:`D_d` is
+the diffusion constant along dimension :math:`d` of the AWH bias.
+For one dimension, :math:`L^2/2D` is the average time to diffuse
+over a distance of :math:`L`. We then takes the maximum crossing
+time over all dimensions involved in the bias.
+Essentially, this formula tells us that a slower system (small :math:`D`)
requires more samples (larger :math:`N^0`) to attain the same level of
accuracy (:math:`\varepsilon_0`) at a given sampling rate. Conversely,
for a system of given diffusion, how to choose the initial biasing rate
maximum free energy difference of the PMF estimate. If this is much
larger than the expected magnitude of the free energy barriers that
should be crossed, then the system is probably being pulled too hard and
-:math:`D` should be decreased. :math:`\varepsilon_0` on the other hand,
-would only be tweaked when starting an AWH simulation using a fairly
-accurate guess of the PMF as input.
+:math:`D` should be decreased. An accurate estimate of the diffusion
+can be obtaining from an AWH simulation with the :ref:`gmx awh` tool.
+:math:`\varepsilon_0` on the other hand, should be a rough estimate
+of the initial error.
Tips for efficient sampling
^^^^^^^^^^^^^^^^^^^^^^^^^^^
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff