[gmx-users] range of magnitude for force constant

[chris.neale at utoronto.ca]

I see two potential problems that might arise with very large force constants.

1. If you start far from equilibrium your system will most definitely  
crash. It appears that this is not going to be an issue for you here  
though.

2. The less mass in your groups and the higher the force constant, the  
faster the oscillations will become. This may force you to reduce your  
timestep in order to obtain a stable system. C-C force constants are  
on the order of 3E+05 kJ/mol/nm^2, and are stable with 5fs timesteps  
(at least) so I would not be surprised at all if a 1E+06 kj/mol/nm^2  
force constant works fine for you if your groups are larger than one  
heavy atom. If you are doing position restraints on individual atoms  
then you might need to utilize a 1fs timestep if you go as high as  
3E+05 kJ/mol/nm^2 due to the hydrogens. This is all just guesswork  
though, I do not have any numerical results for you.

Such high force constants have no buisness giving you enough overlap  
for doing wham with any reasonable number of separate umbrellas. The  
only reason this is going to work for you is that these bilayers are a  
little messy and it is difficult to determine exactly where the  
bilayer reference position is. Therefore, the fact that your protein  
is 2-3A up or down from the center of your restraint does of course  
reduce the resolution of your plot, it should not affect the general  
magnitudes of the free energy profile and increasing the force  
constant may not increase the resolution as much as you expect since  
the noise of the bilayer itself not being homogeneous on Z.


Unless you have access to enormous cpu resources, I suggest staying  
with 500-1000kj/mol/nm^2 and running for longer time at this reduced  
number of windows. Even peptides take a long time to equilibrate in a  
bilayer environment.

Chris.

--- original message ---

I have rather elementary questions regarding force constants used in
restraining position/orientation:

    - I have imposed a force with constants  kx, ky, and kz =10^6
    kJ/mol/nm^2  on a particular atom of a peptide to harmonically position
    restrain the molecule inside a bilayer. Is it alright to use such high
    values for force constants? I have used kz =10^3 kJ/mol/nm^2 but it seems
    that the peptide can still move 2-3 angstrom up and down during the course
    of 200ps simulation so I needed to restrain it more tightly.
    - If I calculate PMF along the z-coordinate (performing separate MD's for
    the peptide along z axis) -- assuming i get reasonable overlaps for the
    histograms -- would such a large number create (numerical) problems in WHAM
    calculations?
    - what is a typical force constant to use in angle_restraint_z (referring
    to formula 4.68 in 3.3 manual). Considering the units of k_c, is 500kJ/mol
    way off?