[gmx-users] position restraint and umbrella potential contribution to the virial

[Christopher Neale]

Dear users and developers:

We have conducted some NPT umbrella sampling (US) simulations to investigate the free energy associate with the aproach of 2 alpha helical peptides. The order parameter is the Cartesian X dimension. During these simulations, in order to maintain (i) the structure of isolated helices and (ii) the relative orientation of the isolated helices, we also applied position restraints in the Cartesian Y and Z dimensions. In addition to free energies of interaction, we are also interested in looking at the system volume as a function of separation distance. In doing so, we realized that the restraint potentials are affecting the virial and thus the volume of the simulation system. I consider this volume dependence to be an artefact of the procedure that we used and we would therefore like to debias this volume data. 

Debiasing is relatively straightforward for the umbrella potential. When dx is less than dx0, the umbrella potential provides an expansion force that increases the volume. Conversely, when dx is greater than dx0, the volume is decreased. However, there is no simple way to debias the volume for the effects of position restraints, which are unpredictable a priori as the helices come into contact. It should be possible to go back to the .xtc files to compute the restraint forces and use these, together with the compressibility, to compute the volume change that was induced by the forces via the virial.

My first question is why these restraint forces are added to the virial at all. Perhaps the US potential needs to be added (although I am not entirely sure why), but if the system is translationally invariant then why would position restraints be included in the virial?

My second question is what kind of errors might I introduce by modifying the source code such that these restraint forces are not included in the virial?

My third question is what gromacs source code file contains the virial contribution of the position restraints? (the forces from the US potential are added to the virial in mdlib/pull.c , where I can set pull->bVirial=FALSE for do_pull_pot() ).

Thank you,
Chris.