[gmx-users] Rotational PMF using angle_restraints_z and virtual definition

[chris.neale at utoronto.ca]

Hello,

I am attempting to construct the free energy profile for rotation of a  
short peptide hairpin. For this purpose, I am using [  
angle_restraints_z ]. I will focus this post on the simplified test  
system that I have been using to characterize the problem that I am  
experiencing: single peptide in vacuo using the sd integrator. The  
free energy profile for rotation here should be flat (After entropic  
debiasing).

However, I do not find a flat PMF. I find that the PMF (generated by  
WHAM) has local maxima at every center of restraint. This indicates  
over-sampling at distances away from each center of restraint and if I  
input a lower force constant to WHAM (400 kJ/mol in place of the 500  
kJ/mol that I used for sampling), I can indeed get a flat PMF. So it  
seems to me that there is some problem with either i) the sampling or  
with ii) the post-sampling analysis of sampled positions. I have ruled  
out some obvious possibilities for option ii: a) The stored angles  
accurately reproduce the sampled angles as checked by comparison to  
g_energy; b) The modified version of WHAM that I am using  
(angle_restraints_z is non-harmonic) successfully produces a flat PMF  
with fake data generated based on a MC generation program that I  
wrote. Further, the local-maxima-at-centers-of-restraint problem is  
not alleviated by increasing the window overlap.

Importantly, I'm using a highly unusual virtual atom architecture to  
get at this angle to the z-axis because I was unable to define a COM  
group for this purpose. I do end up defining the angle between two  
COMs, but I do it via i) generation of a number of [ virtual_sites2 ] :

[ virtual_sites2 ]
;sites from funct a1
    288    52   239     1   0.5

Where I eventually utilize virtual atom 288 as one atom in the
definition of angle_restraints_z. Note that atoms 52 and 239 are not virtual.

[ angle_restraints_z ]
138 288 1 0 500 1

I have created a test system in which I take a water molecule and  
place a virtual atom between the two hydrogens and do similar [  
angle_restraints_z ] on the angle between the oxygen and the new  
virtual atom -- here I do not see the effect.

I have also tested not using the virtual site in my peptide  
simulation, using instead:

[ angle_restraints_z ]
138 239 1 0 500 1

and in this case I also do not see the effect.

And ideas about what could be going on or test systems that I might try?

My best guess is that by putting the restraint on a virtual_sites2  
atom, I am not doing exactly what I think I am doing, but I have no  
idea what the difference might be.

Thank you,
Chris.