Density of states and thermodynamics from energy distributions at different temperatures

This is a tutorial for a set of MATLAB scripts which takes as input potential energy values observed for different replicas in a temperature Replica Exchange (REX) simulation and does the following:

• Computes histograms of potential energy values from each replica
• Plots energy histograms
• Uses energy distributions to compute the partition function at each REX temperature and the density of states based on the Multi Histogram method (details below).
• Uses density of states to compute thermodynamic quantities (Free Energy, Mean Energy, Entropy, Heat Capacity) as a function of temperature
• Plots thermodynamic quantities as a function of temperature

Sample input files and required MATLAB scripts -

ener2thermD.m, DOS2thermod.m and multihist_simple.m

are in clust:/sharedscratch/ss2029/ener2thermoD.

Sample output generated using publishing feature of MATLAB using data from a REX simulation of Alanine Dipeptide (vacuum, AMBER FF) is on my website at http://people.pwf.cam.ac.uk/ss2029/html/ener2thermoD.html

The REX simulations were done using Gromacs 4.5.4 on CLUST. Contact me for details and scripts of the REX simulation on CLUST.

Note on multihist_simple.m:

This implementation of the multi histogram method is based on Section 4A and appendix of Weerasinghe and Amar JCP (1993) [1]. The starting point of the method is a set of linear equations obtained by taking the derivative of $\chi ^{2}$ (= weighted sum of residue squared) as per standard least squares fitting procedure.

Derivative wrt density of state for each bin (Eqn 29 of W&A) gives nB (=number of energy bins) equations:

$V_{k}[\sum _{j}p(k,j)]-\sum _{j}p(k,j)Y_{k}=\sum _{j}p(k,j)W(k,j)\;,k=1,..,nB$

and derivative wrt partition function gives nT (=number of temperatures) equations:

$\sum _{i}[p(i,k)V_{i}]-Y_{k}=\sum _{i}p(i,k)W(i,k)\;,k=1,..,nT$,

where

$V_{i}=\ln[\Omega _{i}dE]\;,i=1,..,nB$, and

$Y_{i}=\ln[Q_{i}]\;,i=1,..,nT$ are unknowns, and

$p(b,t)$ = (normalized) probability of observing energy in bin b at temperature t (specified by input matrix enerPdfs).

W&A first solve for the partition functions by substituting density of state from the first set of equations into second, and then solve for density of states. This implementation sets up the nbins+ntemps linear equations and solves them directly in one go.

Early version of this code was based on Mark Miller's fortran code which implemented the W&A approach.

--Ss2029 13:01, 23 May 2011 (UTC)