From: JC Gumbart (gumbart_at_ks.uiuc.edu)
Date: Thu Mar 01 2012 - 12:49:09 CST
Dear Aron,
The statistical mechanics of rigid waters is different from that of flexible waters. The energy in the Maxwell-Boltzmann distribution for a given water molecule is no longer sum(1/2*mi*vi^2) for all three atoms, but instead 1/2*m_total*v_oxygen^2 (similar for other rigid hydrogens and their parent atoms). The velocities reported for the hydrogens by NAMD, I believe, are meaningless. When you use this formulation, you recover the proper distribution and temperature for the system and its components. I did find a slight difference in the temperature of the solute and water (the solute was about 5% lower), but I don't believe it is significant. The water has a lot more particles, and therefore generates more converged statistical averages (the tutorial even states that agreement only within +/- 10 K is expected).
Also note that the units have not changed in the velocity output by NAMD at any time. I suspect you only obtained the proper temperature for the full system because you empirically derived a conversion factor that made the data fit what you expected.
You might also look at the code that NAMD itself uses to calculate temperature, based on the equipartition theorem, noting in particular that it removes constrained degrees of freedom.
It is certainly worthwhile to verify for yourself the underpinnings of an MD code, but be careful about stating there is something wrong prematurely!
JC
On Feb 28, 2012, at 9:58 PM, Aron Broom wrote:
> 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 : Tue Dec 31 2013 - 23:21:42 CST