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 curvefitting tool however and therefore GROMACS only supports a limited number of functions. Table 18 lists the available options with the corresponding commandline options. The underlying routines for fitting use the LevenbergMarquardt algorithm as implemented in the lmfit package 162 (a barebones version of which is included in GROMACS in which an option for errorweighted 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}}+(1a_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 timecorrelated 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, SteenSæ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{1a_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 stretchedexponential approximation, \(\tau_s\) and \(\beta_s\) are constants for slow relaxation and \(C\) is the prefactor 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.