Analysis

Temperature

Compare the distribution of the instantaneous temperature estimator

\[\mathcal{T}(t) = \frac{2}{3N}\sum_{i=1}^{N}\frac{1}{2}m_i \mathbf{v}_{i}(t)^{2}\]

to the theoretical distribution in the canonical ensemble,

\[P(\mathcal{T}) \propto \exp\left[-\frac{3N(\mathcal{T} - T)^2}{2T^2}\right]\]

Notes on how to do this:

• extract the time series \(\mathcal{T}(t)\) from the energy (EDR) file with g_energy
• histogram the time series with the numpy.histogram() function (if you normalize the distribution with the density=True keyword then the next step will be easier...)
• plot the exact distribution together with the data (you probably want to normalize the exact distribution so that its integral is 1: calculate the normalization constant)
• try to extract \(T\) from your simulation data by fitting the exact distribution. You can use SciPy’s scipy.optimize.leastsq() function. A starting point is shown below.

Least square fitting in Python

SciPy’s scipy.optimize.leastsq() function finds a set of \(M\) parameters \(\{p_k\}_{1\leq k \leq M}\) that minimizes the square of the deviations of a fit function from the \(N\) data points \((x_i, y_i)\). The residual is

\[\delta_i = y_i - f(x_i; \{p_k\})\]

and the function that is minimized over all parameters is the sum of the squared residuals

\[\Delta(\{p_k\}) = \sum_{i=1}^{N} [y_i - f(x_i; \{p_k\})]^2\]

The example shows how to fit a straight line, i.e. the fit function is simply

\[y = f(x) = p_0 x + p_1\]

It has two parameters \(p_0, p_1\). The code only requires you to define the residual in the function errfunc()

\[\delta_i = y_i - (p_0 x_i + p_1)\]

You can write your own function fit_P_NVT() by using the following code as a template

# simple linear fitting with SciPy
import scipy.optimize

def fit_linear(x, y, pinitial):
    """Fit a straight line to the data in *x* and *y*.

    The fitted function is ``y = m*x + t``.

    :Arguments:

      *x*
         array of x values
      *y*
         array of y values
      *pinitial*
         initial guess [slope, offset] e.g. ``[1.0, 0.0]``; this
         must always be a list

    :Returns: Parameters ``[m, t]``
    """
    def errfunc(p,x,y):
        return  y - (p[0]*x + p[1])     # residuals

    p, msg = scipy.optimize.leastsq(errfunc, pinitial[:], args=(x, y))
    print msg
    return p

Before programming, think about the following:

• What data do you want to fit?
• What is your fit function?
• How many parameters does the fit function contain?

Your output should be similar to the example below.

Example output

The figure below shows the output from 100 ps of NVT MD with a coupling time of 1 ps. The tempoerature was recorded every 0.5 ps and thus the histogram is rather jagged. Nevertheless, it agrees with the expected canonical distribution \(P(\mathcal{T})\) and a fit of the distribution to the histogrammed data recovers the average temperature over the simulation (309.9 K (fit) vs 310.6±8.9 K).

Distribution of instantaneous temperature from NVT simulation with Nosé-Hoover thermostat (coupling time 1 ps). Simulation cell with 506 TIP4P water molecules.

Conserved Nosé-Hoover energy

Extract the conserved quantity from the energy (EDR) file with g_energy and plot a time series. Is it conserved?

Additional observables

If you have time, also analyze the following:

• timeseries for

  • total energy E
  • temperature T
  • conserved quantities, if there are any (should all be in the energy (EDR) file but possibly also look in the log file, md.log)

  Use g_energy (and possibly tools like grep or Python to extract data from log).

• square of the fluctuations (i.e. the variance <(X-<X>)**2> for

  • kinetic energy
  • potential energy

• histograms

  • total energy
  • conserved quantities (possibly look in the log file, md.log)

• self diffusion coefficient (g_msd)

• radial pair distribution function between the oxygens (OW-OW)

  To get the OW create an index file with make_ndx:

  make_ndx -f md.tpr -o OW.ndx <<EOF
  keep 0
  del 0
  a OW
  q
  EOF

  Then run g_rdf with only on the selected OW atoms:

  g_rdf -n OW.ndx -s md.tpr -f md.xtc  -xvg none -o rdf.xvg

  Plot rdf.xvg in matplotlib.

  Integrate \(4\pi r g(r) dr\) to the first minimum: how many waters are in the first hydration shell?