.. -*- encoding: utf-8 -*-

Analysis
========

Temperature
-----------

Compare the distribution of the instantaneous temperature estimator

.. math::

   \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,

.. math::

   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 :math:`\mathcal{T}(t)` from the energy (EDR)
  file with :program:`g_energy`
* histogram the time series with the :func:`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 :math:`T` from your simulation data by fitting the
  exact distribution. You can use SciPy's
  :func:`scipy.optimize.leastsq` function. A starting point is shown below.


 

Least square fitting in Python
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

.. math::

   \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

.. math::

   \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

.. math::

   y = f(x) = p_0 x + p_1

It has two parameters :math:`p_0, p_1`. The code only requires you to
define the residual in the function :func:`errfunc`

.. math::

   \delta_i = y_i - (p_0 x_i + p_1)

You can write your own function :func:`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 :math:`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).

.. figure:: /figures/temperature.*
   :width: 30%
   :alt: Distribution of instantaneous temperature from NVT simulation

   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
:program:`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 :program:`g_energy` (and possibly tools like :program:`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 (:program:`g_msd`)

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

  To get the OW create an index file with :program:`make_ndx`::

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

  Then run :program:`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 :math:`4\pi r g(r) dr` to the first minimum: how many
  waters are in the first hydration shell?



.. def P_NVT(T, N, T0=310):
..    sigma2 = T**2/(3*N)
..    return 1/np.sqrt(2*np.pi*sigma2) * np.exp(-(T-T0)**2/(2*sigma2))
..
.. def fit_P_NVT(x, y, N, pinitial):
..    def errfunc(p,x,y, N=N):
..        return  P_NVT(x, N, T0=p[0]) - y
..
..    p, msg = scipy.optimize.leastsq(errfunc, pinitial[:], args=(x, y))
..    print msg
..    return p