[gmx-users] Restraints influeincing free energy calculations
Hi all, Just a quick question - and perhaps some wisdom from those with more experience. I am trying to use metadynamics with gromacs to plot an energy landscape about the rotation of two proteins. If I were to define my collective variables as ThetaA and ThetaB, Theta being the rotation about the protein principle axis, I could sample identical rotations but at different distances apart. Therefore, this introduces the idea of a third collective variable, absolute distance between the centre-of-mass of the proteins. Given that I don't actually want to sample space where these two proteins are beyond peptide-peptide interactions, it is of little interest plus would be extremely computationally crippling, I thought I could apply a restraint between the two proteins at a distance previously identified in literature using FRET. My two questions: 1) Would applying a restraint mess with or bias free energy calculations in any form of free energy calculations, be it, metadynamics or umbrella sampling, and, 2) words of wisdom from anyone who has undertaken such a problem? Many thanks Anthony -- gmx-users mailing listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Frozen covalent bound atoms
Hi all, I would imagine this has been covered before, yet I don't think I have unearthed the right search inquiry yet. I want to make a dihedral angle along the length of a helical protein using the vectors (A, B, C), (B, C, D), where B is a C-alpha on the helix backbone, and A , C and D are fake atoms. I need the fake atom A to be covalently bound/attached to C-alpha atom B, maintaining the constant angle A-B-C. I then need to keep a totally frozen position of fake atoms C and D, which are covalently bound in the order B-C-D. All fake atoms must not interact with any atoms at all. I don't want someone to tell me how to do it, I would just rather someone told me a key word/link/section in the gromacs manual. Many thanks Anthony -- gmx-users mailing listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
RE: [gmx-users] renumbering residue numbers for version 4.0.2
I figured it might of been something like that so I took out the second protein but I was still getting the same error message. Really odd. I can give examples of the numbering, if that helps: Topology order: helixA 1 helixB 1 POPC395 SOL 23765 CL 10 Gro file: All atoms according to the topology file atoms numbered 1 to n. Residue numbers are 1 to k. In order of residue number helixA is 1 to 33 (of which I have made to reflect in the helixA.itp file), and 34 to 67 for helixB (as reflected in helixB.itp) Start of helixA.itp ; nr type resnr residue atom cgnr charge mass typeB chargeB massB ; residue 1 ARG rtp ARG q +1.0 1 NL1ARG N 1 -0.6614.0067 ; qtot -0.66 Start of helixB.itp ; nr type resnr residue atom cgnr charge mass typeB chargeB massB ; residue 34 ARG rtp ARG q +1.0 1 NL34ARG N 1 -0.6614.0067 ; qtot -0.66 I'm thinking of taking everything out of helixB.itp and putting it in helixA.itp, to make something like dimer.itp. But then that would involve making massive changes to the atom ids. Thanks Anthony From: gmx-users-boun...@gromacs.org [gmx-users-boun...@gromacs.org] on behalf of Mark Abraham [mark.j.abra...@gmail.com] Sent: 29 March 2013 22:20 To: Discussion list for GROMACS users Subject: Re: [gmx-users] renumbering residue numbers for version 4.0.2 Looks like something is having trouble dealing with the second protein segment. I would not expect the .itp residue numbering to be significant, except that it changes between residues. Maybe the .gro and .itp residue numbering had to match, back in that prehistoric era... Without seeing what numbers things have, it's hard to say anything useful. Mark On Fri, Mar 29, 2013 at 6:26 PM, Nash, Anthony wrote: > Thanks for that, I've renumbered the residues in the .gro file, and even > so far as in the respective protein .itp files. However, grompp is still > throwing out the same error. Not sure where it is holding onto the number > 652 (which is the start amino acid index of my protein) as I have reset the > residue numbers to begin from 1. > > Residue numbers in the .top are not numbered consecutively from 1 (rather, > resnumber = 652 and resnr_diff = -1) > > My top file is: > > #include "forcefield_popc.itp" > #include "mut_helixA_renumber.itp" > #include "mut_helixB_renumber.itp" > > #include "popc.itp" > #ifdef FLEX_SPC > #include "flexspc.itp" > #else > #include "spc.itp" > #endif > #include "ions.itp" > > [ system ] > Motif A oncogene in water > > [ molecules ] > helixA 1 > helixB 1 > POPC395 > SOL 23765 > CL 10 > > The grompp command is: > ~/gromacs-localpressure/bin/grompp_d -c starting_config.gro -n lateral.ndx > -p ../../mut_system_renumber.top -o local_pressure -f local_pressure.mdp > > Any thoughts? > Many thanks > > > From: gmx-users-boun...@gromacs.org [gmx-users-boun...@gromacs.org] on > behalf of Mark Abraham [mark.j.abra...@gmail.com] > Sent: 29 March 2013 14:17 > To: Discussion list for GROMACS users > Subject: Re: [gmx-users] renumbering residue numbers for version 4.0.2 > > On Fri, Mar 29, 2013 at 2:10 PM, Nash, Anthony > wrote: > > > Apologies if this has been covered once before, I've had a hearty search > > and I've yet to find anything. > > > > I have just installed a modified verison of Gromacs 4.0.2 which calculate > > lateral pressure. However, when I run grommp I get: > > > > Residue numbers in the .top are not numbered consecutively from 1 > (rather, > > resnumber = 652 and resnr_diff = -1) > > > > Yet, running grompp on the same files with version 4.5.1 everything is > > perfectly fine. I am guessing somewhere in between the two versions > > something was changed. > > > > So, does this sound about right? Does the order of residue numbers no > > longer mater with the latest versions of grompp? And ultimately, does > > anyone know how I can go through a .gro file and renumber every residue > > (short of writing a perl script)? > > > > Use genconf -renumber from probably any version (and curse whoever thought > that was a good place to put it :-)) > > Mark > -- > gmx-users mailing listgmx-users@gromacs.org > 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
RE: [gmx-users] renumbering residue numbers for version 4.0.2
Thanks for that, I've renumbered the residues in the .gro file, and even so far as in the respective protein .itp files. However, grompp is still throwing out the same error. Not sure where it is holding onto the number 652 (which is the start amino acid index of my protein) as I have reset the residue numbers to begin from 1. Residue numbers in the .top are not numbered consecutively from 1 (rather, resnumber = 652 and resnr_diff = -1) My top file is: #include "forcefield_popc.itp" #include "mut_helixA_renumber.itp" #include "mut_helixB_renumber.itp" #include "popc.itp" #ifdef FLEX_SPC #include "flexspc.itp" #else #include "spc.itp" #endif #include "ions.itp" [ system ] Motif A oncogene in water [ molecules ] helixA 1 helixB 1 POPC395 SOL 23765 CL 10 The grompp command is: ~/gromacs-localpressure/bin/grompp_d -c starting_config.gro -n lateral.ndx -p ../../mut_system_renumber.top -o local_pressure -f local_pressure.mdp Any thoughts? Many thanks From: gmx-users-boun...@gromacs.org [gmx-users-boun...@gromacs.org] on behalf of Mark Abraham [mark.j.abra...@gmail.com] Sent: 29 March 2013 14:17 To: Discussion list for GROMACS users Subject: Re: [gmx-users] renumbering residue numbers for version 4.0.2 On Fri, Mar 29, 2013 at 2:10 PM, Nash, Anthony wrote: > Apologies if this has been covered once before, I've had a hearty search > and I've yet to find anything. > > I have just installed a modified verison of Gromacs 4.0.2 which calculate > lateral pressure. However, when I run grommp I get: > > Residue numbers in the .top are not numbered consecutively from 1 (rather, > resnumber = 652 and resnr_diff = -1) > > Yet, running grompp on the same files with version 4.5.1 everything is > perfectly fine. I am guessing somewhere in between the two versions > something was changed. > > So, does this sound about right? Does the order of residue numbers no > longer mater with the latest versions of grompp? And ultimately, does > anyone know how I can go through a .gro file and renumber every residue > (short of writing a perl script)? > Use genconf -renumber from probably any version (and curse whoever thought that was a good place to put it :-)) Mark -- gmx-users mailing listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] renumbering residue numbers for version 4.0.2
Apologies if this has been covered once before, I've had a hearty search and I've yet to find anything. I have just installed a modified verison of Gromacs 4.0.2 which calculate lateral pressure. However, when I run grommp I get: Residue numbers in the .top are not numbered consecutively from 1 (rather, resnumber = 652 and resnr_diff = -1) Yet, running grompp on the same files with version 4.5.1 everything is perfectly fine. I am guessing somewhere in between the two versions something was changed. So, does this sound about right? Does the order of residue numbers no longer mater with the latest versions of grompp? And ultimately, does anyone know how I can go through a .gro file and renumber every residue (short of writing a perl script)? Many thanks Anthony -- gmx-users mailing listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
RE: [gmx-users] g_wham missing option
Hi Justin, Thanks for the reply. You were spot on about the version difference. Thanks again. Anthony From: gmx-users-boun...@gromacs.org [gmx-users-boun...@gromacs.org] on behalf of Justin Lemkul [jalem...@vt.edu] Sent: 02 February 2013 18:06 To: Discussion list for GROMACS users Subject: Re: [gmx-users] g_wham missing option On 2/2/13 12:58 PM, Nash, Anthony wrote: > Hi All, > > I am using Gromacs 4.5.5 and running free energy calculations. According to > the article "g_wham - A Free Weighted Histogram.." (J. Chem. Theory. > Comput. 2010, 6, 3713-3720), there is the option -bs-method. However, I am > unable to find this option when running g_wham. I want to use bayesian > bootstraps of complete histograms. > > I have a feeling I am missing something completely obvious. Any help would be > appreciated. > Take a closer look at g_wham -h. The option you're asking about is definitely there. Also be sure you're using the version you think you are; g_wham was overhauled between 4.0.7 and 4.5. -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 listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] g_wham missing option
Hi All, I am using Gromacs 4.5.5 and running free energy calculations. According to the article "g_wham - A Free Weighted Histogram.." (J. Chem. Theory. Comput. 2010, 6, 3713-3720), there is the option -bs-method. However, I am unable to find this option when running g_wham. I want to use bayesian bootstraps of complete histograms. I have a feeling I am missing something completely obvious. Any help would be appreciated. Thanks Anthony -- gmx-users mailing listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
RE: [gmx-users] PMF Transmembrane proteins
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: gmx-users-boun...@gromacs.org [gmx-users-boun...@gromacs.org] on behalf of Justin Lemkul [jalem...@vt.edu] 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
RE: [gmx-users] PMF Transmembrane proteins
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. 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. 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!? 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. 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. I'll get around to putting up some ASAP. Thanks very much for your help Anthony From: gmx-users-boun...@gromacs.org [gmx-users-boun...@gromacs.org] on behalf of Justin Lemkul [jalem...@vt.edu] Sent: 30 December 2012 19:25 To: Discussion list for GROMACS users Subject: Re: [gmx-users] PMF Transmembrane proteins On 12/30/12 12:28 PM, Nash, Anthony wrote: > Dear gmx users, > > I posted a couple of weeks ago with regards to correctly using umbrella > sampling and the WHAM on atomistic transmembrane proteins with a reaction > coordinate as a function of interhelical distance. I have a single TM dimer, > but with a different transmembrane domain face at the helix-helix interface. > So you have proteins with identical sequences in each case, but different orientations? > What did I conclude from the replies? Correct calculation of the differences > in free energy is taken by normalizing to zero at the plateau (flat) of the > graph i.e., where your PMF graph (after g_wham) flattens out you apply the > -zprof0 at this point, before calculating the difference between of each > system. > > The problem I face is that this occurs around 6.8 - 7.5nm along the reaction > coordinate. And, applying a umbrella window every half angstrom has made this > extremely computationally expensive. However, I have persisted, and I have Every half Angstrom? That sounds like massive overkill for almost any system. Do you mean half nanometer instead? > all my graphs, and I have a plateau on all of them. Yet the PMF graphs > between 4 - 7.5 nm are not yet converging or close to the same sampled > energy, even though all four systems are identical in amino acid composition. > Each window has ran for up to 40 ns. > Are there any differences at all in these systems, or are they simply intended to be replicates of one another? I'm just curious to know what you're comparing here. > So when I normalise, depending on where I normalise I will get massively > different difference in free energies. > Real numbers would help us understand what's going on here, as would a report of the error estimate that g_wham cna provide. > I have tried the strategy of taking the windows of a particular systems at > 4.5 nm to 7.5 nm along the reaction coordinate (where there is no short range > interaction), and applying those windows into the other three systems. This > brings their region of plateau a lot closer together, whilst preserving the > windows and the co
[gmx-users] PMF Transmembrane proteins
Dear gmx users, I posted a couple of weeks ago with regards to correctly using umbrella sampling and the WHAM on atomistic transmembrane proteins with a reaction coordinate as a function of interhelical distance. I have a single TM dimer, but with a different transmembrane domain face at the helix-helix interface. What did I conclude from the replies? Correct calculation of the differences in free energy is taken by normalizing to zero at the plateau (flat) of the graph i.e., where your PMF graph (after g_wham) flattens out you apply the -zprof0 at this point, before calculating the difference between of each system. The problem I face is that this occurs around 6.8 - 7.5nm along the reaction coordinate. And, applying a umbrella window every half angstrom has made this extremely computationally expensive. However, I have persisted, and I have all my graphs, and I have a plateau on all of them. Yet the PMF graphs between 4 - 7.5 nm are not yet converging or close to the same sampled energy, even though all four systems are identical in amino acid composition. Each window has ran for up to 40 ns. So when I normalise, depending on where I normalise I will get massively different difference in free energies. I have tried the strategy of taking the windows of a particular systems at 4.5 nm to 7.5 nm along the reaction coordinate (where there is no short range interaction), and applying those windows into the other three systems. This brings their region of plateau a lot closer together, whilst preserving the windows and the convergence of the region along the reaction coordinate where there are close-range forces in play (converges 0.8 nm - 4.5 nm). Alternatively, could I make the assumption that as their amino acid composition is identical across all four systems, I can normalise where each of the four graphs finish converging (around 4.5 nm)? I would really appreciate any advice from someone who has insight on this issue. I have looked through many journals and founds lots of entries where they just stop at 2.0 nm with no plateau, yet I haven't found a single WHAM of reconstructing umbrella windows where they normalise at the plateau. I have also contacted a few authors of journals with regards to why they cut off at 2 nm. None of them took into consideration that their PMF graphs were still raising. Many thanks Anthony -- gmx-users mailing listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
RE: [gmx-users] Where to stop with TM protein PMF calculations
Hi Justin, Thanks for the swift reply. I am now convinced that the 'plateau' approach is the best. This brings me to a further concern. Take two systems (umbrella sampling along a reaction coordinate etc..) which are identical in composition, yet different in conformation e.g., a helical tilt/crossing angle. I run the two systems along a reaction coordinate until I have a plateau. Both systems exhibit a plateau at around 6.5 - 7.5 nm. However, there is a significant difference in free energy (as a theoretical example 50 kj mol^-1 difference) at this region, which would bias the difference in free energy when I normalise both graphs to zero at this point. Given that both systems are atomistically identical, one would assume an identical free energy when the helices are far enough apart that they don't interact. By using g_wham in increments of 1 ns (0-1 ns, 0-2 ns, 0-3 ns,.), I can see how the curve near the interfacing helices converges, yet at 4+ nm the curve is still dropping. Would you advice that I run additional windows near the plateau region so it is identical in value across both systems? Many thanks Anthony From: gmx-users-boun...@gromacs.org [gmx-users-boun...@gromacs.org] on behalf of Justin Lemkul [jalem...@vt.edu] Sent: 19 December 2012 13:23 To: Discussion list for GROMACS users Subject: Re: [gmx-users] Where to stop with TM protein PMF calculations On 12/19/12 4:12 AM, Nash, Anthony wrote: > Good morning, > > A bit of a long one I am afraid. > > I am simulating a transmembrane dimer, and calculating the association free > energy through potential of mean force calculations as a function of > interhelical distance. I have got very good umbrella coverage along my > reaction coordinate, however, I would like to know where I should stop > calculating and normalise to zero? I am comparing two systems of identical > composition but with a different conformation. Hence, the normalising step is > vital. > > Looking in the literature the cut offs used seem to be around 2nm (two JACS > papers come to mind). I've pulled as far as needed to reach a plateau in the > PMF curve. This is around 6.5 nm to 8 nm. The problem is if I normalise at > the plateau (which seems to be the beliefs of the professors at my > institute), the comparative free energy between the two systems is very > different to normalising at 2nm (note: the papers in literature do not > observe a plateau, there is still an obvious upward trend at their cutoff). > > I need justification of where to normalise. According to literature they cut > off outside of interaction range. I have decided to test umbrella windows > along the reaction coordinate by decomposing the energy of the system. LJ > short and Coul short between peptides at a reaction coordinate distance of 8 > nm, is 0 as expected. However, I am still getting energy values on Coul. > recip even when I decompose (although the gromacs literature says I can't) > the system by systematically setting every charge to zero, i.e., > > > Peptide A (A), Peptide B (B), SOL, lipid and counter ions have zero charge > throughout. > E_coul_recip_A_B = E_coul_recip_(A_B,A_A,B_B) - (E_coul_recip_A_A + > E_coul_recip_B_B) > > All charges are zero, with the exception to peptide A in E_coul_recip_A_A and > peptide B in E_coul_recip_B_B. and peptide A and B in > E_coul_recip_(A_B,A_A,B_B). When I add E_coul_recip_A_B to the LJ short and > Coul short I do not get an energy of zero. > > So: > 1) The coul-recip decomposition, is this a flawed approach? The gromacs > manual says it can't be decomposed. It's not trivial to do. There are some posts about really detailed ways that you might pull it off, but I don't know whether it's worth the effort. If the umbrella sampling simulations were conducted with PME, there is never a true "non-interacting" state. There is always some finite contribution to long-range electrostatics terms. > 2) Where along the reaction coordinate should I cut off? I would believe the longer distance. At 2 nm, you still may have induced order in either lipids or water between your two proteins, thus affecting the forces between them. I see no reason to think that if a plateau has not been established that the results should be correct. It could be argued that all you're doing is picking some random point along the reaction coordinate and calling it "dissociated" for the sake of convenience. > 3) Why don't authors pull until they reach a plateau? > I have seen several published papers that appear to have very arbitrary stopping points, so just because it's published, doesn't mean it's perfect :) I believe in a plateau, a demonstration that, within the cavea
[gmx-users] Where to stop with TM protein PMF calculations
Good morning, A bit of a long one I am afraid. I am simulating a transmembrane dimer, and calculating the association free energy through potential of mean force calculations as a function of interhelical distance. I have got very good umbrella coverage along my reaction coordinate, however, I would like to know where I should stop calculating and normalise to zero? I am comparing two systems of identical composition but with a different conformation. Hence, the normalising step is vital. Looking in the literature the cut offs used seem to be around 2nm (two JACS papers come to mind). I've pulled as far as needed to reach a plateau in the PMF curve. This is around 6.5 nm to 8 nm. The problem is if I normalise at the plateau (which seems to be the beliefs of the professors at my institute), the comparative free energy between the two systems is very different to normalising at 2nm (note: the papers in literature do not observe a plateau, there is still an obvious upward trend at their cutoff). I need justification of where to normalise. According to literature they cut off outside of interaction range. I have decided to test umbrella windows along the reaction coordinate by decomposing the energy of the system. LJ short and Coul short between peptides at a reaction coordinate distance of 8 nm, is 0 as expected. However, I am still getting energy values on Coul. recip even when I decompose (although the gromacs literature says I can't) the system by systematically setting every charge to zero, i.e., Peptide A (A), Peptide B (B), SOL, lipid and counter ions have zero charge throughout. E_coul_recip_A_B = E_coul_recip_(A_B,A_A,B_B) - (E_coul_recip_A_A + E_coul_recip_B_B) All charges are zero, with the exception to peptide A in E_coul_recip_A_A and peptide B in E_coul_recip_B_B. and peptide A and B in E_coul_recip_(A_B,A_A,B_B). When I add E_coul_recip_A_B to the LJ short and Coul short I do not get an energy of zero. So: 1) The coul-recip decomposition, is this a flawed approach? The gromacs manual says it can't be decomposed. 2) Where along the reaction coordinate should I cut off? 3) Why don't authors pull until they reach a plateau? Many many thanks Anthony -- gmx-users mailing listgmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
