docs/reference-manual/analysis/curve-fitting.rst

   1 Curve fitting in |Gromacs|
   2 --------------------------
   3
   4 Sum of exponential functions
   5 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   6
   7 Sometimes it is useful to fit a curve to an analytical function, for
   8 example in the case of autocorrelation functions with noisy tails.
   9 |Gromacs| is not a general purpose curve-fitting tool however and
  10 therefore |Gromacs| only supports a limited number of functions.
  11 :numref:`Table %s <table-fitfn>`  lists the available options with the corresponding
  12 command-line options. The underlying routines for fitting use the
  13 Levenberg-Marquardt algorithm as implemented in the lmfit package \ :ref:`162 <reflmfit>`
  14 (a bare-bones version of which is included in |Gromacs| in which an
  15 option for error-weighted fitting was implemented).
  16
  17 .. |exp|  replace:: :math:`e^{-t/{a_0}}`
  18 .. |aexp| replace:: :math:`a_1e^{-t/{a_0}}`
  19 .. |exp2| replace:: :math:`a_1e^{-t/{a_0}}+(1-a_1)e^{-t/{a_2}}`
  20 .. |exp5| replace:: :math:`a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_4`
  21 .. |exp7| replace:: :math:`a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_5e^{-t/{a_4}}+a_6`
  22 .. |exp9| replace:: :math:`a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_5e^{-t/{a_4}}+a_7e^{-t/{a_6}}+a_8`
  23 .. |nexp2| replace:: :math:`a_2\ge a_0\ge 0`
  24 .. |nexp5| replace:: :math:`a_2\ge a_0\ge 0`
  25 .. |nexp7| replace:: :math:`a_4\ge a_2\ge a_0 \ge0`
  26 .. |nexp9| replace:: :math:`a_6\ge a_4\ge a_2\ge a_0\ge 0`
  27
  28 .. _table-fitfn:
  29
  30 .. table:: Overview of fitting functions supported in (most) analysis tools
  31     that compute autocorrelation functions. The **Note** column describes
  32     properties of the output parameters.
  33     :align: center
  34     :widths: auto
  35
  36     +-------------+------------------------------+---------------------+
  37     | Command     | Functional form :math:`f(t)` | Note                |
  38     | line option |                              |                     |
  39     +=============+==============================+=====================+
  40     | exp         | |exp|                        |                     |
  41     +-------------+------------------------------+---------------------+
  42     | aexp        | |aexp|                       |                     |
  43     +-------------+------------------------------+---------------------+
  44     | exp_exp     | |exp2|                       | |nexp2|             |
  45     +-------------+------------------------------+---------------------+
  46     | exp5        | |exp5|                       | |nexp5|             |
  47     +-------------+------------------------------+---------------------+
  48     | exp7        | |exp7|                       | |nexp7|             |
  49     +-------------+------------------------------+---------------------+
  50     | exp9        | |exp9|                       | |nexp9|             |
  51     +-------------+------------------------------+---------------------+
  52
  53
  54 Error estimation
  55 ~~~~~~~~~~~~~~~~
  56
  57 Under the hood |Gromacs| implements some more fitting functions, namely a
  58 function to estimate the error in time-correlated data due to Hess \ :ref:`149 <refHess2002a>`:
  59
  60 .. math:: \varepsilon^2(t) =
  61           2 \alpha\tau_1\left(1+\frac{\tau_1}{t}\left(e^{-t/\tau_1}-1\right)\right)
  62                 + 2 (1-\alpha)\tau_2\left(1+\frac{\tau_2}{t}\left(e^{-t/\tau_2}-1\right)\right)
  63           :label: eqntimecorrerror
  64
  65 where :math:`\tau_1` and :math:`\tau_2` are time constants (with
  66 :math:`\tau_2 \ge \tau_1`) and :math:`\alpha` usually is close to 1 (in
  67 the fitting procedure it is enforced that :math:`0\leq\alpha\leq 1`).
  68 This is used in :ref:`gmx analyze <gmx analyze>` for error estimation using
  69
  70 .. math:: \lim_{t\rightarrow\infty}\varepsilon(t) = \sigma\sqrt{\frac{2(\alpha\tau_1+(1-\alpha)\tau_2)}{T}}
  71           :label: eqnanalyzeerrorest
  72
  73 where :math:`\sigma` is the standard deviation of the data set and
  74 :math:`T` is the total simulation time \ :ref:`149 <refHess2002a>`.
  75
  76 Interphase boundary demarcation
  77 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  78
  79 In order to determine the position and width of an interface,
  80 Steen-Sæthre *et al.* fitted a density profile to the following function
  81
  82 .. math:: f(x) ~=~ \frac{a_0+a_1}{2} - \frac{a_0-a_1}{2}{\rm
  83           erf}\left(\frac{x-a_2}{a_3^2}\right)
  84           :label: eqndesprofilefunc
  85
  86 where :math:`a_0` and :math:`a_1` are densities of different phases,
  87 :math:`x` is the coordinate normal to the interface, :math:`a_2` is the
  88 position of the interface and :math:`a_3` is the width of the
  89 interface \ :ref:`163 <refSteen-Saethre2014a>`. This is implemented
  90 in :ref:`gmx densorder <gmx densorder>`.
  91
  92 Transverse current autocorrelation function
  93 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  94
  95 In order to establish the transverse current autocorrelation function
  96 (useful for computing viscosity  \ :ref:`164 <refPalmer1994a>`) the following function is
  97 fitted:
  98
  99 .. math:: f(x) ~=~ e^{-\nu}\left({\rm cosh}(\omega\nu)+\frac{{\rm
 100           sinh}(\omega\nu)}{\omega}\right)
 101           :label: eqntransverseautocorrfunc
 102
 103 with :math:`\nu = x/(2a_0)` and :math:`\omega = \sqrt{1-a_1}`. This is
 104 implemented in :ref:`gmx tcaf <gmx tcaf>`.
 105
 106 Viscosity estimation from pressure autocorrelation function
 107 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 108
 109 The viscosity is a notoriously difficult property to extract from
 110 simulations \ :ref:`149 <refHess2002a>`, :ref:`165 <refWensink2003a>`. It is *in principle*
 111 possible to determine it by integrating the pressure autocorrelation
 112 function \ :ref:`160 <refPSmith93c>`, however this is often hampered by
 113 the noisy tail of the ACF. A workaround to this is fitting the ACF to
 114 the following function \ :ref:`166 <refGuo2002b>`:
 115
 116 .. math:: f(t)/f(0) = (1-C) {\rm cos}(\omega t) e^{-(t/\tau_f)^{\beta_f}} + C
 117           e^{-(t/\tau_s)^{\beta_s}}
 118           :label: eqnviscestpressureautocorr
 119
 120 where :math:`\omega` is the frequency of rapid pressure oscillations
 121 (mainly due to bonded forces in molecular simulations), :math:`\tau_f`
 122 and :math:`\beta_f` are the time constant and exponent of fast
 123 relaxation in a stretched-exponential approximation, :math:`\tau_s` and
 124 :math:`\beta_s` are constants for slow relaxation and :math:`C` is the
 125 pre-factor that determines the weight between fast and slow relaxation.
 126 After a fit, the integral of the function :math:`f(t)` is used to
 127 compute the viscosity:
 128
 129 .. math:: \eta = \frac{V}{k_B T}\int_0^{\infty} f(t) dt
 130           :label: eqncompviscosity
 131
 132 This equation has been applied to computing the bulk and shear
 133 viscosity using different elements from the pressure tensor \ :ref:`167 <refFanourgakis2012a>`.