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

===========
 Solutions
===========


Nosé-Hoover thermostat in MDP file
==================================

Your MDP file should contain the following lines in order to use the
Nosé-Hoover thermostat::
   
   ; TEMPERATURE COUPLING
   Tcoupl                   = nose-hoover
   tc-grps                  = SYSTEM
   tau-t                    = 1
   ref-t                    = 310

It's also a good idea to start with initial velocities at the same
temperature::

   gen-temp                 = 310


Canonical distribution of the temperature estimator
===================================================

The normalized distribution of the temperature estimator is 

.. math::

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

Extract the temperature as a function of time from the energy file ::
 
  echo "Temperature" | g_energy -f md.edr -o energy.xvg -xvg none

and histogram the temperature time series . The code below shows how
to do this and also includes 

* the function :func:`P_NVT` to plot the normalized canonical
  distribution :math:`P(\mathcal{T})`

* and the code to fit the normalized histogram to
  :math:`P(\mathcal{T})` for comparison

One can run it and it will read the data from :file:`energy.xvg` and
produce the plots :file:`temperature.{pdf|png|svg}`. (In the code,
:math:`\mathcal{T}` is written ``T`` and :math:`T` is ``T0``.)::

   #!/usr/bin/env python
   # -*- encoding: utf-8 -*-
 
   import numpy as np
   import scipy.optimize

   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

   if __name__ == "__main__":
       import matplotlib.pyplot as plt

       N_particles = 506

       # load energy file containing time, temperature
       time, T = np.loadtxt("energy.xvg", unpack=True)
       print("simulation timeseries: T+/-sigma(T) = {0:.2f}+/-{1:.2f} K".format(T.mean(), T.std()))

       # normalized histogram from simulation
       h, e = np.histogram(T, bins=20, density=True)
       m = 0.5*(e[1:]+e[:-1])

       # use <T> as starting point for fit
       T_fit, = fit_P_NVT(m, h, N_particles, [T.mean()])
       print("fit to P(T): T_fit = {0:.2f} K".format(T_fit))

       # plotting
       trange = np.linspace(280, 350, 100)

       fig = plt.figure(figsize=(5,5))
       ax = fig.add_subplot(111)

       ax.plot(m, h, 'k-', label="simulation")
       ax.plot(trange, P_NVT(trange, N_particles, T0=310), 'b-', label=r"canonical $P(T)$")
       ax.plot(trange, P_NVT(trange, N_particles, T0=T_fit), 'r--', lw=2, label=r"fit, $T=309.9$ K")
       ax.set_xlabel("temperature $T$ (K)")
       ax.legend(loc="best")
       ax.set_title(u"Nosé-Hoover thermostat at $T=310$ K")
       fig.savefig("temperature.pdf")
       fig.savefig("temperature.png")
       fig.savefig("temperature.svg")