1Curve fitting in |Gromacs| 2-------------------------- 3 4Sum of exponential functions 5~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 6 7Sometimes it is useful to fit a curve to an analytical function, for 8example in the case of autocorrelation functions with noisy tails. 9|Gromacs| is not a general purpose curve-fitting tool however and 10therefore |Gromacs| only supports a limited number of functions. 11:numref:`Table %s <table-fitfn>` lists the available options with the corresponding 12command-line options. The underlying routines for fitting use the 13Levenberg-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 15option 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 54Error estimation 55~~~~~~~~~~~~~~~~ 56 57Under the hood |Gromacs| implements some more fitting functions, namely a 58function 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 65where :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 67the fitting procedure it is enforced that :math:`0\leq\alpha\leq 1`). 68This 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 73where :math:`\sigma` is the standard deviation of the data set and 74:math:`T` is the total simulation time \ :ref:`149 <refHess2002a>`. 75 76Interphase boundary demarcation 77~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 78 79In order to determine the position and width of an interface, 80Steen-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 86where :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 88position of the interface and :math:`a_3` is the width of the 89interface \ :ref:`163 <refSteen-Saethre2014a>`. This is implemented 90in :ref:`gmx densorder <gmx densorder>`. 91 92Transverse current autocorrelation function 93~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 94 95In order to establish the transverse current autocorrelation function 96(useful for computing viscosity \ :ref:`164 <refPalmer1994a>`) the following function is 97fitted: 98 99.. math:: f(x) ~=~ e^{-\nu}\left({\rm cosh}(\omega\nu)+\frac{{\rm 100 sinh}(\omega\nu)}{\omega}\right) 101 :label: eqntransverseautocorrfunc 102 103with :math:`\nu = x/(2a_0)` and :math:`\omega = \sqrt{1-a_1}`. This is 104implemented in :ref:`gmx tcaf <gmx tcaf>`. 105 106Viscosity estimation from pressure autocorrelation function 107~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 108 109The viscosity is a notoriously difficult property to extract from 110simulations \ :ref:`149 <refHess2002a>`, :ref:`165 <refWensink2003a>`. It is *in principle* 111possible to determine it by integrating the pressure autocorrelation 112function \ :ref:`160 <refPSmith93c>`, however this is often hampered by 113the noisy tail of the ACF. A workaround to this is fitting the ACF to 114the 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 120where :math:`\omega` is the frequency of rapid pressure oscillations 121(mainly due to bonded forces in molecular simulations), :math:`\tau_f` 122and :math:`\beta_f` are the time constant and exponent of fast 123relaxation in a stretched-exponential approximation, :math:`\tau_s` and 124:math:`\beta_s` are constants for slow relaxation and :math:`C` is the 125pre-factor that determines the weight between fast and slow relaxation. 126After a fit, the integral of the function :math:`f(t)` is used to 127compute the viscosity: 128 129.. math:: \eta = \frac{V}{k_B T}\int_0^{\infty} f(t) dt 130 :label: eqncompviscosity 131 132This equation has been applied to computing the bulk and shear 133viscosity using different elements from the pressure tensor \ :ref:`167 <refFanourgakis2012a>`. 134