Hi Justin, Thanks for the help. The block analysis has done wonders! I had been interested in seeing how the PMF graphs converged, i.e., 0-1ns, 0-2ns, 0-3ns. Of course this always includes the earliest time steps. Long story short, by cutting off the first 10ns so I am only sampling 10ns-40ns, the free energy curves are very different.
Regions where the peptides are close together are not that different. But between 3nm and 7.5nm the curves plateau with a lot less energy difference between the like-with-like system. I will run bootstrap analysis, get out the error, see how that goes. Oddly enough the curve is not as smooth now, even though there is ample histogram coverage. I may push the simulations for longer. Thanks, I think this is now moving forwards again. Anthony ________________________________________ From: [email protected] [[email protected]] on behalf of Justin Lemkul [[email protected]] Sent: 30 December 2012 20:08 To: Discussion list for GROMACS users Subject: Re: [gmx-users] PMF Transmembrane proteins On 12/30/12 2:57 PM, Nash, Anthony wrote: > Dear Justin > > Thanks for your reply. > > I am studying a TM peptide, looking at how favourable the front and reverse > face are. Hence, different orientations but the same composition of amino > acids. > > I then have a duplicate of this system (call it system B), but with a couple > of residue substitutions. I am comparing free energies of System A front > orientation with System A back orientation (and then B with B). > > I have separated the peptides, and ran umbrella simulations every 1 angstrom > apart (I experimented with windows closer together, in some regions hence the > generalisation to half an angstrom). Papers I have looked through which > sample up to 2 nm apart have been running umbrella windows every 0.2 of an > angstrom. My PMF plateau is around 6.5-7.5 nm, which implies around 78 > windows per system. > Without seeing the actual histograms and .mdp settings, I can't comment on this approach directly. It seems like you should be able to use far fewer windows, perhaps with a slightly less restrictive force constant to allow for slightly broader distributions. That's a hand-waving guess, though, since I've seen none of the data. > Taking one system as an example, if I normalise at 7.5 nm the difference at > their global minima is around 37.63 kJ mol−1 (8.99 kCal mol−1), but then when > I break down the PMF graph (0-1ns, 0-2ns, 0-3nm) the regions greater than 4.5 > nm have not converged. Hence, normalising anywhere in this unconverged region > would be wildly inaccurate. If I normalise to zero at the point in which one > of the system pairs (A&A, or B&B) converges, the free energy different is 5 > kJ mol−1. > I think it would be more productive to do your block analysis such that you discard initial frames to see what effect the (presumably) unequilibrated part of the trajectory is. Even with prior equilibration (i.e. with position restraints), you still need to leave some frames behind as equilibration. Analyzing 0-40 ns, 10-40 ns, 20-40 ns, etc will probably be more informative. If you're always considering all frames in the WHAM analysis, that can throw off the results considerably. > If by error estimate you mean from the inbuilt bootstrap analysis, then I am > very confident that the error along the entire reaction coordinate is very > small; which is confusing as I am not sure whether the slower converging > regions of the reaction coordinate (4.5 - 7.5) should shower greater error!? > That's not uncommon. At greater distance, you have two peptides floating around rather independently, so the restraint potential may vary a bit more. At closer COM distances, the interactions between the peptides keep the orientations a bit more stable. > Yes it is heart breaking to hear that I may need to run the system for > longer. This however, isn't a massive problem, given the time of the year the > HPCs I am using are a bit empty. Would you suggest I need to run the 40 ns > for longer (continuation run), or begin brand new windows from scratch and > include those in the g_wham calculations? I guess the 40 ns one, given that > lipid reorientation needs to be taken into account. > There is no need to scrap the existing data. You may only need to extend the simulations, but do carry out the analysis I suggested before to see if there is undue influence from the initial frames. > I don't have numbers in front of me, but the system box size is big (I went > through your tutorial several times, and I have seen the PMF graph as a > result of a box too small). The longest i.e., the reaction coordinate, is > along the Y axes, Z is normal to the bilayer. When I calculate the PMF, I > have pull_dim = Y Y N. > Good to know. Lipid ordering effects can extend up to a few nm though, so that's why I asked. -Justin -- ======================================== Justin A. Lemkul, Ph.D. Research Scientist Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ======================================== -- gmx-users mailing list [email protected] http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to [email protected]. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- gmx-users mailing list [email protected] http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to [email protected]. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
