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).

Table 18 Overview of fitting functions supported in (most) analysis tools that compute autocorrelation functions. The Note column describes properties of the output parameters.#

Command line option

Functional form f(t)

Note

exp

et/a0

aexp

a1et/a0

exp_exp

a1et/a0+(1a1)et/a2

a2a00

exp5

a1et/a0+a3et/a2+a4

a2a00

exp7

a1et/a0+a3et/a2+a5et/a4+a6

a4a2a00

exp9

a1et/a0+a3et/a2+a5et/a4+a7et/a6+a8

a6a4a2a00

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:

(450)#ε2(t)=2ατ1(1+τ1t(et/τ11))+2(1α)τ2(1+τ2t(et/τ21))

where τ1 and τ2 are time constants (with τ2τ1) and α usually is close to 1 (in the fitting procedure it is enforced that 0α1). This is used in gmx analyze for error estimation using

(451)#limtε(t)=σ2(ατ1+(1α)τ2)T

where σ 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

(452)#f(x) = a0+a12a0a12erf(xa2a32)

where a0 and a1 are densities of different phases, x is the coordinate normal to the interface, a2 is the position of the interface and a3 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:

(453)#f(x) = eν(cosh(ων)+sinh(ων)ω)

with ν=x/(2a0) and ω=1a1. 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:

(454)#f(t)/f(0)=(1C)cos(ωt)e(t/τf)βf+Ce(t/τs)βs

where ω is the frequency of rapid pressure oscillations (mainly due to bonded forces in molecular simulations), τf and βf are the time constant and exponent of fast relaxation in a stretched-exponential approximation, τs and β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:

(455)#η=VkBT0f(t)dt

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