# Estimate Partition Function from a simulation

I have a system and I've carried out a long molecular dynamics simulation over it. I would like to estimate the partition function $Z.$ Theoretically, one would compute: $$Z=\dfrac{1}{N!h^{3N}}\int \exp \left(-\beta\frac{-p^2}{2m} \right) \exp\left(-\beta V(r) \right) \ dp \ dr,$$ but if the system I'm dealing with is, say, a protein in water or something like that, this integral is of course absolutely untractable. When thinking of how to estimate the partition function of susch a system from a molecular dynamics simulation, I guess the first thing that comes to mind is the following: Consider all the snapshots obtained during your simulation, label them with the index $j$ and call $E_j$ to the energy of each snapshot. For giant thermalized systems, I guess we can neglect the probability of a concrete configuration repeating, so I believe a nice try would be: $$\tilde{Z}=\frac{dV_{\text{ps}}}{N!h^{3N}}\sum_{j} \exp\left( -\beta E_j \right),$$ where $dV_{\text{pc}}$ is the volume of small cube in phase space (its volume would have to be adjusted manually depending oh how many snapshots one has taken).

My question, briefly put, is, how can I estimate the error of estimating $Z$ by $\tilde{Z}$?

It seems intuitive that typically $\tilde{Z}$ should be a good aproximation to $Z,$ but one can easilly think of situations where this isn't true: Imagine that my initial configuration lies in, so as to speak, a valley in the energetic landscape, surrounded by a chain of "high mountains" (high-energy configuration) and that at the other side of these "mountains" lie low-energy configuration. It is clear that in such a case there could be a considerable error when estimating $Z$ from $\tilde{Z}.$ But since one does not know the aspect of the energetic landscape beforehand, it seems difficult to deal with this.

The question stops here, but just to provide context I'll explain briefly what I'm thinking about: I want to study a chemical reaction, and to that end I would like to know the free energy profile along a certain reaction coordinate. For that, an Umbrella Sampling technique is used: you drag the reaction coordinate along the reaction path, forcing it to be at certain positions at different steps (otherwhise, the simulation wouldn't take you over the energetic barriers in short simulation times). I am aware of several methods to compute free energies, like FEP or WHAM, which is quite sophisticated, but I was guessing that if the technique I mentioned above to estimate the partition function $Z$ was correct, then the goal of getting the free energy profile would be fairly simple: you just drag along the reaction coordinate, perform molecular dynamics simulations for different restricted values of the reaction coordinate (you force the reaction coordinate to be at a concrete position by means of a harmonic potential), and at each value of the reaction coordinate you estimate the partition function as above. That's the idea and that's why I am interested in having some idea on methods to estimate the error of approximating $Z$ by $\tilde{Z}.$

Note: this question was originally in Physics.stackexchange. It didn't receive any attention at all and I asked for it to be migrated to Chemistry.stackexchange, but the flag was ignored, so I've decided to delete my question on Physics.stacexchange and post it here

## 2 Answers

What you're saying does make sense. You can, in fact directly estimate the partition function using nested sampling, with error bounds. See for example this paper http://journals.aps.org/prb/abstract/10.1103/PhysRevB.93.174108 . In your case without PBC it's actually even simpler. Also check out Livia's blog here https://liviabpartay.wordpress.com/publications/

There's a deeper problem, in my opinion: your quantity \bar{Z} will not go to the partition function in the limit of infinite sampling. The problem is that your configurations j are already drawn from a Boltzmann distribution (if your simulation is at a defined T), and your expression for \bar{Z} then weights these Boltzmann-distributed energies by a second Boltzmann factor.

Getting the partition function (or configuration integral) from a thermalized simulation is not easy, and hence is not typically attempted.

You might be interested in a discussion of this in Andrew Leach's book, Eqs 6.21-6.28, in the 2001 edition.

• Yes, this is precisely what I was worried about. I understand you. Weighting the energies with the Boltzmann factors would make sense if the different snapshots were uniform in the configuration space, but the sampling is already biased towars more probable regions. So, to estimate the partition function, I would need the density of states... Is the Wand and Landau algorithm the best way to do that? – user28429 Sep 15 '16 at 13:54