# Curve fitting in GROMACS¶

## Sum of exponential functions¶

Sometimes it is useful to fit a curve to an analytical function, for example in the case of autocorrelation functions with noisy tails. GROMACS is not a general purpose curve-fitting tool however and therefore GROMACS only supports a limited number of functions. Table 18 lists the available options with the corresponding command-line options. The underlying routines for fitting use the Levenberg-Marquardt algorithm as implemented in the lmfit package 162 (a bare-bones version of which is included in GROMACS in which an option for error-weighted fitting was implemented).

Command line option | Functional form \(f(t)\) | Note |
---|---|---|

exp | \(e^{-t/{a_0}}\) | |

aexp | \(a_1e^{-t/{a_0}}\) | |

exp_exp | \(a_1e^{-t/{a_0}}+(1-a_1)e^{-t/{a_2}}\) | \(a_2\ge a_0\ge 0\) |

exp5 | \(a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_4\) | \(a_2\ge a_0\ge 0\) |

exp7 | \(a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_5e^{-t/{a_4}}+a_6\) | \(a_4\ge a_2\ge a_0 \ge0\) |

exp9 | \(a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_5e^{-t/{a_4}}+a_7e^{-t/{a_6}}+a_8\) | \(a_6\ge a_4\ge a_2\ge a_0\ge 0\) |

## Error estimation¶

Under the hood GROMACS implements some more fitting functions, namely a function to estimate the error in time-correlated data due to Hess 149:

where \(\tau_1\) and \(\tau_2\) are time constants (with \(\tau_2 \ge \tau_1\)) and \(\alpha\) usually is close to 1 (in the fitting procedure it is enforced that \(0\leq\alpha\leq 1\)). This is used in gmx analyze for error estimation using

where \(\sigma\) is the standard deviation of the data set and \(T\) is the total simulation time 149.

## Interphase boundary demarcation¶

In order to determine the position and width of an interface,
Steen-Sæthre *et al.* fitted a density profile to the following function

where \(a_0\) and \(a_1\) are densities of different phases, \(x\) is the coordinate normal to the interface, \(a_2\) is the position of the interface and \(a_3\) is the width of the interface 163. This is implemented in gmx densorder.

## Transverse current autocorrelation function¶

In order to establish the transverse current autocorrelation function (useful for computing viscosity 164) the following function is fitted:

with \(\nu = x/(2a_0)\) and \(\omega = \sqrt{1-a_1}\). This is implemented in gmx tcaf.

## Viscosity estimation from pressure autocorrelation function¶

The viscosity is a notoriously difficult property to extract from
simulations 149, 165. It is *in principle*
possible to determine it by integrating the pressure autocorrelation
function 160, however this is often hampered by
the noisy tail of the ACF. A workaround to this is fitting the ACF to
the following function 166:

where \(\omega\) is the frequency of rapid pressure oscillations (mainly due to bonded forces in molecular simulations), \(\tau_f\) and \(\beta_f\) are the time constant and exponent of fast relaxation in a stretched-exponential approximation, \(\tau_s\) and \(\beta_s\) are constants for slow relaxation and \(C\) is the pre-factor that determines the weight between fast and slow relaxation. After a fit, the integral of the function \(f(t)\) is used to compute the viscosity:

This equation has been applied to computing the bulk and shear viscosity using different elements from the pressure tensor 167.