rigidBonds water and langevinHydrogen off are INCOMPATIBLE when using a solute with hdyrogens

From: Aron Broom (broomsday_at_gmail.com)
Date: Tue Feb 28 2012 - 21:58:25 CST

Dear NAMD users,

I've looked around on the mailing list and couldn't find a clear answer to
as to whether or not it is supposed to be ok to use rigidBonds water (say
with TIP3P) but then leave the solute (say a protein) flexible with a 1 fs
timestep and langevinHydrogen off in order to save on calculating
collisions with the large number of rigid hydrogens in the sample. It
might seem like one should not use langevinHydrogen off as long as there
are any flexible hydrogens, but several sources, including an NAMD tutorial
(
http://www.ks.uiuc.edu/Training/Tutorials/science/channel/channel-tutorial-files/ABF/win1/win1A.conf)
have exactly this setup.

Well, I believe I have confirmed that you should NOT do this. In looking
at a 60ns simulation of ~2000 protein atoms in an ~100,000 water atom box,
and determining the temperature of the non-rigid parts of the system by
calculating their kinetic energies (KE=0.5mv^2), binning them, and then
fitting to the boltzmann distribution ( y = (2/sqrt(Pi * (KbT)^3 )) *
sqrt(x) * exp(-x/KbT) ), I find that the temperature of the whole system is
fine (set to 298 K, and comes out on average at 297 K), but the temperature
of the protein alone is ~20-25 K higher than the surrounding solvent
throughout the simulation.

This is a big problem, especially for studying ligand binding, and I
imagine for anything having a 20-25 K difference across the solvent-solute
interface is bad, not to mention just being incorrect.

I tested running an extra 4ns from the end of my 60 ns, with two changes in
NAMD parameters:

1) just turn langevinHydrogen on, which comes with a 5% performance
penalty, and

2) leave langevinHydrogen off, but change rigidBonds water to rigidBonds
all, and change from a 1 fs timestep to a 2 fs timestep, in which you gain
~30% performance, but some people are sceptical of the results.

Both of these changes alone are capable of bringing the system back to a
reasonable configuration in which the protein and solvent have
approximately the same temperatures (the protein was still slightly hotter
(~5 K, but perhaps a gamma of 1/ps is not fast enough to see complete
equilibrium in 4ns).

Has anyone else seen this and can confirm that the setup of using
ridigBonds water and langevinHydrogen off are not to be used together when
you have a solute that has hydrogens, despite this being somewhat common
practice? If anyone wants to test the temperatures of various parts of
their systems from existing *.vel files, I just followed the instructions
on this page (
http://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/node13.html),
but note that the units in the .vel files must changed since this was
written, so you'll have to convert amu to kg, and angstrom/ps to m/s, and
then use the correct boltzmann constant for those units.

~Aron

-- 
Aron Broom M.Sc
PhD Student
Department of Chemistry
University of Waterloo

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:21:15 CST