Could someone elaborate on how GROMACS’s test-particle insertion subroutine calculates the values found in tpi.xvg and tpidist.xvg? More specifically, how should the energies from tpidist.xvg relate to e.g. the excess chemical potential in an NpT ensemble found in the first column of tpi.xvg? For instance, by taking the values X = beta U - ln(V/) from tpidist.xvg and calculate thereupon <exp(-X)> = <exp(-beta U + ln(V/))>_N does not yield the value exp( - mu_ex / R T ). Unfortunately I cannot find anything in the documentation regarding the content of tpidist.xvg and tpi.xvg.
Thanks for your reply. I have seen the legend of the tpi.xvg file and it is quite self-explanatory indeed. Maybe my initial question was a bit vague. My question is regarding the relation between the first column in tpi.xvg (<mu_ex> in an NpT ensemble) and the data from tpidist.xvg:
Is it true that one can calculate the excess chemical potential \mu_\mathrm{ex} from the tpidist.xvg data by evaluating
\langle \mu_\mathrm{ex} \rangle = R T \ln \langle \exp[-(\beta U -\ln (V/\langle V \rangle)] \rangle_\mathrm{N},
where the quantity in the exponential, the insertion energy \beta U - \ln V/\langle V \rangle, is the data from the first column in tpidist.xvg. The averaging is then obtained by using the absolute counts as weights. Here I assumed that we only consider one single frame with N insertions. I have included \log(V/\langle V \rangle) because GROMACS seems to do that, but it should be equal to zero because we only consider one frame (and thus \langle V \rangle = V, independent of the ensemble).
The first column contains the average over all frames up to the current one. So the value in first column for the last frame is equal to the finally computed chemical potential.
Nit: mu_ex should not have angular brackets as it’s an ensemble property by definition. TPI computes an approximation of mu_ex.
Is the way I would calculate my \mu_\mathrm{ex} from the tpidist.xvg data correct then? Because I don’t get the excess chemical potential that is listed in the tpi.xvg file using the insertion energy distribution. I would just like to double check and make sure that I understand everything.
Thanks, although it seems quite straightforward, I don’t manage to get the value of \mu_\mathrm{ex} from the insertion energy in tpidist.xvg. Just to be 100% sure, I should be able to use the insertion energy as is right? Since
\mu_\mathrm{ex} = -RT \ln \langle \exp\{-\underbrace{(\beta U - \log V/\langle V \rangle)}_\text{energy from tpidist.xvg} \} \rangle = -RT \ln \frac{\langle \exp(-\beta U) V \rangle}{\langle V \rangle}.
Then, I should be able to obtain \langle \exp(-\beta U + \ln V / \langle V \rangle) \rangle by a weighted average of the first column which is the absolute count c_i for the i\mathrm{th} insertion energy U_i^\mathrm{TPI} = \beta U_i - \ln V_i / \langle V \rangle in tpidist.xvg. So in terms of pseudocode something like:
bU = get_bins_of_histogram('tpidist.xvg') # The bins represent: beta U + ln(V/<V>)
counts = get_counts_of_histogram('tpidist.xvg') # get the absolute counts per bin
avg_exp_insert_energy = average( exp(-bU), weights = counts)
mu_ex = -R*T*log( avg_exp_insert_energy ) # J/mol
This should give the correct results (for TPI on a single frame) or do I miss something crucial here?
I see now that I actually wrote the correct formulas for the output in the tpidist.xvg, but that my text was incorrect. There is no per frame output in tpidist.xvg, so you can not recompute mu_ex. Only the last line should match mu_ex.