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

\[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 P_NVT() to plot the normalized canonical distribution \(P(\mathcal{T})\)
• and the code to fit the normalized histogram to \(P(\mathcal{T})\) for comparison

One can run it and it will read the data from energy.xvg and produce the plots temperature.pdf|png|svg. (In the code, \(\mathcal{T}\) is written T and \(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")