------ at 5 ps ----- TAUP density E CPU(sec) .2 (default) .6973 -2841 600 .1 .8389 -2865 618 .04 .9604 -2951 670 .02 .9714 -2980 689 .01 .9722 -2957 698 .008 .9635 -2964 703 .001 .6165 -1821 576The files for this comparison are in the equil_notes directory and are named are md0_2.out, md0_1.out, etc. respectively; TAUP=.01 is the canonical md0.out in the main directory.
Looking at the density over time for some of these:
We see that when the coupling is the same as the step length (.001 ps), the density fluctuates wildly. The default coupling of .2 ps converges very slowly. Results for the more interesting TAUP's over 10 ps:
All three are about the same, for practical purposes, and further convergence toward 1.0 is very gradual.
Finally, running the TAUP=.01 trajectory for another 5 ps using TAUP=.1 to see if the density smooths as in TAUP=.1 in the first graph - and while we're at it, let's compare to 15 ps of TAUP=.1 to see if that would have converged:
The apparent smoothness of TAUP=.1 in the first graph seems to be mostly due to the steep rise in density.
My opinion is that in terms of density, at least, 5 ps of TAUP=.01 would suffice for perturbation of the water box. For this small system where 15 ps takes about half an hour on an HP workstation, the effort expended in deciding what to do next far exceeds the cost of just running, say, 50 ps with the default density. However, these results may be useful for larger systems for which 50 ps may take days or weeks. Let's check some other equilibration parameters to make sure the faster density protocol is not hurting some other aspect of the equilibration.
The faster protocol seems about as good as the less-tightly-coupled version, given that when the pressure coupling is relaxed at 10 ps, the behavior more or less matches the looser version. (A common question is "why does pressure fluctuate so much?" - this is answered in one of our Questions pages.)
The faster protocol seems better: at 5 ps, the potential energy has already equilibrated, and relaxing TAUP at 10 ps has no effect.
We have looked at some macroscopic values that are mostly determined by the water bath - but what about the methane? Since SHAKE is used on covalent bonds involving hydrogen (NTC=NTF=2, i.e. all bonds in the system), the only internal reaction that the methane can have to its environment is in angle energy; since angle potential is not evaluated for the waters (they are held rigid by an artificial H-H bond and the use of SHAKE), it happens that the angle energy in the .out files is all methane:
Here it appears that TAUP=.1 for 15 ps may have a slight edge over the faster protocol for 5 ps, since the longer protocol seems to offer lower and less variable energy. Testing whether this makes a difference in the free energy perturbation of the methane is left as an exercise.