Re: [gmx-users] deltaG from PMF
I run analysis of my PMF curve based on the simulation time: http://speedy.sh/6GyDN/PMF-time.jpg I would surely neglect first 20 ns as the drifts are really high. But then which of those curve would you take as the final PMF? As increasing the time the PMF get lower deltaG values, would you take the last 10 ns then? please, help, Steven On Fri, Aug 24, 2012 at 4:13 PM, Steven Neumann s.neuman...@gmail.com wrote: Thanks for this! Steven On Fri, Aug 24, 2012 at 3:26 PM, Thomas Schlesier schl...@uni-mainz.de wrote: As a first step, i would shift all curves so, that the energy of the minium is for all plots at the same (aribarity) value. The minimum should be the point which has sampled the best. If you shift then all values, it should be easier to spot differences between the plots. And probably make the analysis for every 10ns slices 10-20, 20-30, ... then it's easier to see if you have a drift or fluctuations. If you identify problematic regions you could do there additional umbrella simulations, probably with higher force constants, in order to sample that region better. In theory the PMF should converge if you wait long enough, till the system equilibrates under the external restraints (how long this takes? nobody knows). On could probably wait a long time, big question here is what is the biggest error you would tolerate. On the other hand one can't simulate forever... Look for what other people used for simulation-times for umbrella sampling (i have the impression that 50ns is rather long, but i could fool me here) and what they did to estimate if the calculations are converged. Then decide, what to do. Think nobody here would / will tell you this or this error is ok ... but if you do something reasonable, which is in accord what others do / did it should be ok. Am 24.08.2012 14:59, schrieb gmx-users-requ...@gromacs.org: I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25,schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US
Re: [gmx-users] deltaG from PMF
I would simulate longer. What are the autocorrelation times for the underlying data? Erik 28 aug 2012 kl. 16.20 skrev Steven Neumann: I run analysis of my PMF curve based on the simulation time: http://speedy.sh/6GyDN/PMF-time.jpg I would surely neglect first 20 ns as the drifts are really high. But then which of those curve would you take as the final PMF? As increasing the time the PMF get lower deltaG values, would you take the last 10 ns then? please, help, Steven On Fri, Aug 24, 2012 at 4:13 PM, Steven Neumann s.neuman...@gmail.com wrote: Thanks for this! Steven On Fri, Aug 24, 2012 at 3:26 PM, Thomas Schlesier schl...@uni-mainz.de wrote: As a first step, i would shift all curves so, that the energy of the minium is for all plots at the same (aribarity) value. The minimum should be the point which has sampled the best. If you shift then all values, it should be easier to spot differences between the plots. And probably make the analysis for every 10ns slices 10-20, 20-30, ... then it's easier to see if you have a drift or fluctuations. If you identify problematic regions you could do there additional umbrella simulations, probably with higher force constants, in order to sample that region better. In theory the PMF should converge if you wait long enough, till the system equilibrates under the external restraints (how long this takes? nobody knows). On could probably wait a long time, big question here is what is the biggest error you would tolerate. On the other hand one can't simulate forever... Look for what other people used for simulation-times for umbrella sampling (i have the impression that 50ns is rather long, but i could fool me here) and what they did to estimate if the calculations are converged. Then decide, what to do. Think nobody here would / will tell you this or this error is ok ... but if you do something reasonable, which is in accord what others do / did it should be ok. Am 24.08.2012 14:59, schrieb gmx-users-requ...@gromacs.org: I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25,schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin
Re: [gmx-users] deltaG from PMF
50 ns for each window takes around a week on 12 cpus. I think that neglecting first 20 ns would be the best. The penalty is app 0.5 kcal/mol at this time so really good. On Tue, Aug 28, 2012 at 3:28 PM, Erik Marklund er...@xray.bmc.uu.se wrote: I would simulate longer. What are the autocorrelation times for the underlying data? Erik 28 aug 2012 kl. 16.20 skrev Steven Neumann: I run analysis of my PMF curve based on the simulation time: http://speedy.sh/6GyDN/PMF-time.jpg I would surely neglect first 20 ns as the drifts are really high. But then which of those curve would you take as the final PMF? As increasing the time the PMF get lower deltaG values, would you take the last 10 ns then? please, help, Steven On Fri, Aug 24, 2012 at 4:13 PM, Steven Neumann s.neuman...@gmail.com wrote: Thanks for this! Steven On Fri, Aug 24, 2012 at 3:26 PM, Thomas Schlesier schl...@uni-mainz.de wrote: As a first step, i would shift all curves so, that the energy of the minium is for all plots at the same (aribarity) value. The minimum should be the point which has sampled the best. If you shift then all values, it should be easier to spot differences between the plots. And probably make the analysis for every 10ns slices 10-20, 20-30, ... then it's easier to see if you have a drift or fluctuations. If you identify problematic regions you could do there additional umbrella simulations, probably with higher force constants, in order to sample that region better. In theory the PMF should converge if you wait long enough, till the system equilibrates under the external restraints (how long this takes? nobody knows). On could probably wait a long time, big question here is what is the biggest error you would tolerate. On the other hand one can't simulate forever... Look for what other people used for simulation-times for umbrella sampling (i have the impression that 50ns is rather long, but i could fool me here) and what they did to estimate if the calculations are converged. Then decide, what to do. Think nobody here would / will tell you this or this error is ok ... but if you do something reasonable, which is in accord what others do / did it should be ok. Am 24.08.2012 14:59, schrieb gmx-users-requ...@gromacs.org: I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25,schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM,
Re: [gmx-users] deltaG from PMF
Indeed, extended simulations are always painful. Nevertheless, g_wham calculates the autocorrelations and use them for producing more accurate error estimates. They can also give a hint about wether you simulated long enough. Erik 28 aug 2012 kl. 16.42 skrev Steven Neumann: 50 ns for each window takes around a week on 12 cpus. I think that neglecting first 20 ns would be the best. The penalty is app 0.5 kcal/mol at this time so really good. On Tue, Aug 28, 2012 at 3:28 PM, Erik Marklund er...@xray.bmc.uu.se wrote: I would simulate longer. What are the autocorrelation times for the underlying data? Erik 28 aug 2012 kl. 16.20 skrev Steven Neumann: I run analysis of my PMF curve based on the simulation time: http://speedy.sh/6GyDN/PMF-time.jpg I would surely neglect first 20 ns as the drifts are really high. But then which of those curve would you take as the final PMF? As increasing the time the PMF get lower deltaG values, would you take the last 10 ns then? please, help, Steven On Fri, Aug 24, 2012 at 4:13 PM, Steven Neumann s.neuman...@gmail.com wrote: Thanks for this! Steven On Fri, Aug 24, 2012 at 3:26 PM, Thomas Schlesier schl...@uni-mainz.de wrote: As a first step, i would shift all curves so, that the energy of the minium is for all plots at the same (aribarity) value. The minimum should be the point which has sampled the best. If you shift then all values, it should be easier to spot differences between the plots. And probably make the analysis for every 10ns slices 10-20, 20-30, ... then it's easier to see if you have a drift or fluctuations. If you identify problematic regions you could do there additional umbrella simulations, probably with higher force constants, in order to sample that region better. In theory the PMF should converge if you wait long enough, till the system equilibrates under the external restraints (how long this takes? nobody knows). On could probably wait a long time, big question here is what is the biggest error you would tolerate. On the other hand one can't simulate forever... Look for what other people used for simulation-times for umbrella sampling (i have the impression that 50ns is rather long, but i could fool me here) and what they did to estimate if the calculations are converged. Then decide, what to do. Think nobody here would / will tell you this or this error is ok ... but if you do something reasonable, which is in accord what others do / did it should be ok. Am 24.08.2012 14:59, schrieb gmx-users-requ...@gromacs.org: I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or
Re: [gmx-users] deltaG from PMF
Thanks. Could you give the example of such command? e.g. g_wham -ac -nBootstrap 200 -bins 200 ? On Tue, Aug 28, 2012 at 4:19 PM, Erik Marklund er...@xray.bmc.uu.se wrote: Indeed, extended simulations are always painful. Nevertheless, g_wham calculates the autocorrelations and use them for producing more accurate error estimates. They can also give a hint about wether you simulated long enough. Erik 28 aug 2012 kl. 16.42 skrev Steven Neumann: 50 ns for each window takes around a week on 12 cpus. I think that neglecting first 20 ns would be the best. The penalty is app 0.5 kcal/mol at this time so really good. On Tue, Aug 28, 2012 at 3:28 PM, Erik Marklund er...@xray.bmc.uu.se wrote: I would simulate longer. What are the autocorrelation times for the underlying data? Erik 28 aug 2012 kl. 16.20 skrev Steven Neumann: I run analysis of my PMF curve based on the simulation time: http://speedy.sh/6GyDN/PMF-time.jpg I would surely neglect first 20 ns as the drifts are really high. But then which of those curve would you take as the final PMF? As increasing the time the PMF get lower deltaG values, would you take the last 10 ns then? please, help, Steven On Fri, Aug 24, 2012 at 4:13 PM, Steven Neumann s.neuman...@gmail.com wrote: Thanks for this! Steven On Fri, Aug 24, 2012 at 3:26 PM, Thomas Schlesier schl...@uni-mainz.de wrote: As a first step, i would shift all curves so, that the energy of the minium is for all plots at the same (aribarity) value. The minimum should be the point which has sampled the best. If you shift then all values, it should be easier to spot differences between the plots. And probably make the analysis for every 10ns slices 10-20, 20-30, ... then it's easier to see if you have a drift or fluctuations. If you identify problematic regions you could do there additional umbrella simulations, probably with higher force constants, in order to sample that region better. In theory the PMF should converge if you wait long enough, till the system equilibrates under the external restraints (how long this takes? nobody knows). On could probably wait a long time, big question here is what is the biggest error you would tolerate. On the other hand one can't simulate forever... Look for what other people used for simulation-times for umbrella sampling (i have the impression that 50ns is rather long, but i could fool me here) and what they did to estimate if the calculations are converged. Then decide, what to do. Think nobody here would / will tell you this or this error is ok ... but if you do something reasonable, which is in accord what others do / did it should be ok. Am 24.08.2012 14:59, schrieb gmx-users-requ...@gromacs.org: I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for
Re: [gmx-users] deltaG from PMF
28 aug 2012 kl. 17.26 skrev Steven Neumann: Thanks. Could you give the example of such command? e.g. g_wham -ac -nBootstrap 200 -bins 200 ? Yes. -oiact further sets to what file the correlation times will be written. To save a bit of time, -ac-trestart can perhaps be increased from 1 ps to e.g. 10. On Tue, Aug 28, 2012 at 4:19 PM, Erik Marklund er...@xray.bmc.uu.se wrote: Indeed, extended simulations are always painful. Nevertheless, g_wham calculates the autocorrelations and use them for producing more accurate error estimates. They can also give a hint about wether you simulated long enough. Erik 28 aug 2012 kl. 16.42 skrev Steven Neumann: 50 ns for each window takes around a week on 12 cpus. I think that neglecting first 20 ns would be the best. The penalty is app 0.5 kcal/mol at this time so really good. On Tue, Aug 28, 2012 at 3:28 PM, Erik Marklund er...@xray.bmc.uu.se wrote: I would simulate longer. What are the autocorrelation times for the underlying data? Erik 28 aug 2012 kl. 16.20 skrev Steven Neumann: I run analysis of my PMF curve based on the simulation time: http://speedy.sh/6GyDN/PMF-time.jpg I would surely neglect first 20 ns as the drifts are really high. But then which of those curve would you take as the final PMF? As increasing the time the PMF get lower deltaG values, would you take the last 10 ns then? please, help, Steven On Fri, Aug 24, 2012 at 4:13 PM, Steven Neumann s.neuman...@gmail.com wrote: Thanks for this! Steven On Fri, Aug 24, 2012 at 3:26 PM, Thomas Schlesier schl...@uni-mainz.de wrote: As a first step, i would shift all curves so, that the energy of the minium is for all plots at the same (aribarity) value. The minimum should be the point which has sampled the best. If you shift then all values, it should be easier to spot differences between the plots. And probably make the analysis for every 10ns slices 10-20, 20-30, ... then it's easier to see if you have a drift or fluctuations. If you identify problematic regions you could do there additional umbrella simulations, probably with higher force constants, in order to sample that region better. In theory the PMF should converge if you wait long enough, till the system equilibrates under the external restraints (how long this takes? nobody knows). On could probably wait a long time, big question here is what is the biggest error you would tolerate. On the other hand one can't simulate forever... Look for what other people used for simulation-times for umbrella sampling (i have the impression that 50ns is rather long, but i could fool me here) and what they did to estimate if the calculations are converged. Then decide, what to do. Think nobody here would / will tell you this or this error is ok ... but if you do something reasonable, which is in accord what others do / did it should be ok. Am 24.08.2012 14:59, schrieb gmx-users-requ...@gromacs.org: I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this:
Re: [gmx-users] deltaG from PMF
I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesier schl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schrieb gmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure,
[gmx-users] deltaG from PMF
As a first step, i would shift all curves so, that the energy of the minium is for all plots at the same (aribarity) value. The minimum should be the point which has sampled the best. If you shift then all values, it should be easier to spot differences between the plots. And probably make the analysis for every 10ns slices 10-20, 20-30, ... then it's easier to see if you have a drift or fluctuations. If you identify problematic regions you could do there additional umbrella simulations, probably with higher force constants, in order to sample that region better. In theory the PMF should converge if you wait long enough, till the system equilibrates under the external restraints (how long this takes? nobody knows). On could probably wait a long time, big question here is what is the biggest error you would tolerate. On the other hand one can't simulate forever... Look for what other people used for simulation-times for umbrella sampling (i have the impression that 50ns is rather long, but i could fool me here) and what they did to estimate if the calculations are converged. Then decide, what to do. Think nobody here would / will tell you this or this error is ok ... but if you do something reasonable, which is in accord what others do / did it should be ok. Am 24.08.2012 14:59, schrieb gmx-users-requ...@gromacs.org: I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25,schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ;
Re: [gmx-users] deltaG from PMF
Thanks for this! Steven On Fri, Aug 24, 2012 at 3:26 PM, Thomas Schlesier schl...@uni-mainz.de wrote: As a first step, i would shift all curves so, that the energy of the minium is for all plots at the same (aribarity) value. The minimum should be the point which has sampled the best. If you shift then all values, it should be easier to spot differences between the plots. And probably make the analysis for every 10ns slices 10-20, 20-30, ... then it's easier to see if you have a drift or fluctuations. If you identify problematic regions you could do there additional umbrella simulations, probably with higher force constants, in order to sample that region better. In theory the PMF should converge if you wait long enough, till the system equilibrates under the external restraints (how long this takes? nobody knows). On could probably wait a long time, big question here is what is the biggest error you would tolerate. On the other hand one can't simulate forever... Look for what other people used for simulation-times for umbrella sampling (i have the impression that 50ns is rather long, but i could fool me here) and what they did to estimate if the calculations are converged. Then decide, what to do. Think nobody here would / will tell you this or this error is ok ... but if you do something reasonable, which is in accord what others do / did it should be ok. Am 24.08.2012 14:59, schrieb gmx-users-requ...@gromacs.org: I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25,schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md ; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0
[gmx-users] deltaG from PMF
Dear Gmx Users, I got confused about reading free energy difference from PMF curve. It is the difference between maximum value on PMF curve (plateau - final state) and the starting point corresponding to minima. So e.g. my curve starts at 0 [kcal/mol] going staright to minima of -2 kcal/mol and then increase obtaining plateau of 14 kcal/mol. DeltaG = -16 kcal/mol or -12 kcal/mol ? Please clarify me this, Steven -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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] deltaG from PMF
On 8/21/12 9:36 AM, Steven Neumann wrote: Dear Gmx Users, I got confused about reading free energy difference from PMF curve. It is the difference between maximum value on PMF curve (plateau - final state) and the starting point corresponding to minima. So e.g. my curve starts at 0 [kcal/mol] going staright to minima of -2 kcal/mol and then increase obtaining plateau of 14 kcal/mol. DeltaG = -16 kcal/mol or -12 kcal/mol ? It seems that your system has discovered a more energetically favorable location where the minimum is at -2 kcal/mol. It should not be assumed that whatever arbitrary start point we choose is the energy minimum. In your case, it is not. DeltaG is the difference between final and initial states. Only you know what these states are. -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 * Only plain text messages are allowed! * 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] deltaG from PMF
On Tue, Aug 21, 2012 at 3:07 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 9:36 AM, Steven Neumann wrote: Dear Gmx Users, I got confused about reading free energy difference from PMF curve. It is the difference between maximum value on PMF curve (plateau - final state) and the starting point corresponding to minima. So e.g. my curve starts at 0 [kcal/mol] going staright to minima of -2 kcal/mol and then increase obtaining plateau of 14 kcal/mol. DeltaG = -16 kcal/mol or -12 kcal/mol ? It seems that your system has discovered a more energetically favorable location where the minimum is at -2 kcal/mol. It should not be assumed that whatever arbitrary start point we choose is the energy minimum. In your case, it is not. DeltaG is the difference between final and initial states. Only you know what these states are. -Justin Thank you. Ok - mu curve starts at 0 kcal/mol corresponding to the reaction coordinate of 0.5 nm. going to -2 kcal mol at 0.6 nm then increasing as I said before to plateau of 14 kcal/mol. My first window (shown by grompp) corresponds to 0.8 nm of the reaction coordinate... So the initial state should be at 0.8 nm which corresponds to 0.5 kcal/mol or always take the initial state as initial of the curve? Sorry, got confused. Steven -- 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 * Only plain text messages are allowed! * 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 * Only plain text messages are allowed! * 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] deltaG from PMF
On 8/21/12 10:23 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:07 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 9:36 AM, Steven Neumann wrote: Dear Gmx Users, I got confused about reading free energy difference from PMF curve. It is the difference between maximum value on PMF curve (plateau - final state) and the starting point corresponding to minima. So e.g. my curve starts at 0 [kcal/mol] going staright to minima of -2 kcal/mol and then increase obtaining plateau of 14 kcal/mol. DeltaG = -16 kcal/mol or -12 kcal/mol ? It seems that your system has discovered a more energetically favorable location where the minimum is at -2 kcal/mol. It should not be assumed that whatever arbitrary start point we choose is the energy minimum. In your case, it is not. DeltaG is the difference between final and initial states. Only you know what these states are. -Justin Thank you. Ok - mu curve starts at 0 kcal/mol corresponding to the reaction coordinate of 0.5 nm. going to -2 kcal mol at 0.6 nm then increasing as I said before to plateau of 14 kcal/mol. My first window (shown by grompp) corresponds to 0.8 nm of the reaction coordinate... So the initial state should be at 0.8 nm which corresponds to 0.5 kcal/mol or always take the initial state as initial of the curve? That depends on how you've set up your .mdp files. If you're telling a configuration in which the groups are separated by 0.8 nm that they should be separated by 0.5 nm instead, you need to allow for some time for that to adjust. The restraint distances are what you tell them to be. Your initial configurations should be as close to that COM separation as possible. So, the answer is, I don't have an answer, because you didn't provide an .mdp file. -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 * Only plain text messages are allowed! * 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] deltaG from PMF
Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 500 ; kJ mol^-1 nm^-2 pull_nstxout= 1000; every 2 ps pull_nstfout= 1000 ; every 2 ps I will appreciate your help, Steven On Tue, Aug 21, 2012 at 3:27 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 10:23 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:07 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 9:36 AM, Steven Neumann wrote: Dear Gmx Users, I got confused about reading free energy difference from PMF curve. It is the difference between maximum value on PMF curve (plateau - final state) and the starting point corresponding to minima. So e.g. my curve starts at 0 [kcal/mol] going staright to minima of -2 kcal/mol and then increase obtaining plateau of 14 kcal/mol. DeltaG = -16 kcal/mol or -12 kcal/mol ? It seems that your system has discovered a more energetically favorable location where the minimum is at -2 kcal/mol. It should not be assumed that whatever arbitrary start point we choose is the energy minimum. In your case, it is not. DeltaG is the difference between final and initial states. Only you know what these states are. -Justin Thank you. Ok - mu curve starts at 0 kcal/mol corresponding to the reaction coordinate of 0.5 nm. going to -2 kcal mol at 0.6 nm then increasing as I said before to plateau of 14 kcal/mol. My first window (shown by grompp) corresponds to 0.8 nm of the reaction coordinate... So the initial state should be at 0.8 nm which corresponds to 0.5 kcal/mol or always take the initial state as initial of the curve? That depends on how you've set up your .mdp files. If you're telling a configuration in which the groups are separated by 0.8 nm that they should be separated by 0.5 nm instead, you need to allow for some time for that to adjust. The restraint distances are what you tell them to be. Your initial configurations should be as close to that COM separation as possible. So, the answer is, I don't have an answer, because you didn't provide an .mdp file. -Justin -- Justin A. Lemkul, Ph.D. Research
Re: [gmx-users] deltaG from PMF
On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 500 ; kJ mol^-1 nm^-2 pull_nstxout= 1000; every 2 ps pull_nstfout= 1000 ; every 2 ps Based on these settings you're allowing grompp to set the reference distance to whatever it finds in the coordinate file. It seems clear to me that the sampling indicates what I said before - you have an energy minimum somewhere other than where you started with. What that state corresponds to relative to what you think is going on is for you to decide based on the nature of your system. Whatever is occurring at 0.6 nm of COM separation is of particular interest, since the energy minimum is so distinct. I hope you're doing a thorough analysis of convergence if you're generating velocities at the outset of each run, and removing unequilibrated frames from your analysis. -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 * Only plain text messages are allowed! * 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] deltaG from PMF
On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 500 ; kJ mol^-1 nm^-2 pull_nstxout= 1000; every 2 ps pull_nstfout= 1000 ; every 2 ps Based on these settings you're allowing grompp to set the reference distance to whatever it finds in the coordinate file. It seems clear to me that the sampling indicates what I said before - you have an energy minimum somewhere other than where you started with. What that state corresponds to relative to what you think is going on is for you to decide based on the nature of your system. Whatever is occurring at 0.6 nm of COM separation is of particular interest, since the energy minimum is so distinct. So based on this the deltaG will correspond to -5.22 as the initial state was taken at 0.4 nm corresponding to 0 kcal/mol as the moment corresponding to the minimum is the coordinate from SMD where last hydrogen bond was broken. Would you agree? I hope you're doing a thorough analysis of convergence if you're generating velocities at the outset of each run, and removing unequilibrated frames from your analysis. When I use WHAM I skip first 200 ps of each window as the equilibration. Steven -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 * Only plain text messages are allowed! * 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 list
Re: [gmx-users] deltaG from PMF
On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 500 ; kJ mol^-1 nm^-2 pull_nstxout= 1000; every 2 ps pull_nstfout= 1000 ; every 2 ps Based on these settings you're allowing grompp to set the reference distance to whatever it finds in the coordinate file. It seems clear to me that the sampling indicates what I said before - you have an energy minimum somewhere other than where you started with. What that state corresponds to relative to what you think is going on is for you to decide based on the nature of your system. Whatever is occurring at 0.6 nm of COM separation is of particular interest, since the energy minimum is so distinct. So based on this the deltaG will correspond to -5.22 as the initial state was taken at 0.4 nm corresponding to 0 kcal/mol as the moment corresponding to the minimum is the coordinate from SMD where last hydrogen bond was broken. Would you agree? Based on the very little information I have, no. It would appear that the 0.4 nm separation is in fact some metastable state and the true energy minimum is at 0.6 nm of COM separation. What's going on at that location? I hope you're doing a thorough analysis of convergence if you're generating velocities at the outset of each run, and removing unequilibrated frames from your analysis. When I use WHAM I skip first 200 ps of each window as the equilibration. That seems fairly short, especially given the generation of velocities in conjunction with the Parrinello-Rahman barostat, which can be very temperamental. -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 * Only plain text messages are allowed! * Please
Re: [gmx-users] deltaG from PMF
On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 500 ; kJ mol^-1 nm^-2 pull_nstxout= 1000; every 2 ps pull_nstfout= 1000 ; every 2 ps Based on these settings you're allowing grompp to set the reference distance to whatever it finds in the coordinate file. It seems clear to me that the sampling indicates what I said before - you have an energy minimum somewhere other than where you started with. What that state corresponds to relative to what you think is going on is for you to decide based on the nature of your system. Whatever is occurring at 0.6 nm of COM separation is of particular interest, since the energy minimum is so distinct. So based on this the deltaG will correspond to -5.22 as the initial state was taken at 0.4 nm corresponding to 0 kcal/mol as the moment corresponding to the minimum is the coordinate from SMD where last hydrogen bond was broken. Would you agree? Based on the very little information I have, no. It would appear that the 0.4 nm separation is in fact some metastable state and the true energy minimum is at 0.6 nm of COM separation. What's going on at that location? My mistake. The initial grompp of 1st configuartion (where ligand stakced on keratin surface) corresponds to 0.6 nm where is the minimum. Thus deltaG would be -7.22 kcal/mol. Am I right? Or Shall I take difference between 0 and 5.22 ? I hope you're doing a thorough analysis of convergence if you're generating velocities at the outset of each run, and removing unequilibrated frames from your analysis. When I use WHAM I skip first 200 ps of each window as the equilibration. That seems fairly short, especially given the generation of velocities in conjunction with the Parrinello-Rahman barostat, which can be very temperamental. Would you suggest e.g. skip 1
Re: [gmx-users] deltaG from PMF
On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 500 ; kJ mol^-1 nm^-2 pull_nstxout= 1000; every 2 ps pull_nstfout= 1000 ; every 2 ps Based on these settings you're allowing grompp to set the reference distance to whatever it finds in the coordinate file. It seems clear to me that the sampling indicates what I said before - you have an energy minimum somewhere other than where you started with. What that state corresponds to relative to what you think is going on is for you to decide based on the nature of your system. Whatever is occurring at 0.6 nm of COM separation is of particular interest, since the energy minimum is so distinct. So based on this the deltaG will correspond to -5.22 as the initial state was taken at 0.4 nm corresponding to 0 kcal/mol as the moment corresponding to the minimum is the coordinate from SMD where last hydrogen bond was broken. Would you agree? Based on the very little information I have, no. It would appear that the 0.4 nm separation is in fact some metastable state and the true energy minimum is at 0.6 nm of COM separation. What's going on at that location? My mistake. The initial grompp of 1st configuartion (where ligand stakced on keratin surface) corresponds to 0.6 nm where is the minimum. Thus deltaG would be -7.22 kcal/mol. Am I right? Or Shall I take difference between 0 and 5.22 ? -7.22 kcal/mol sounds much more logical to me. If your first configuration is at the energy minimum, that's not something you ignore. The zero point can be set wherever you like with the g_wham flag -zprof0, so it's really rather arbitrary. The WHAM algorithm simply sets the leftmost window (smallest value along the reaction coordinate) to zero to construct the remainder of the PMF curve. I hope you're doing a thorough analysis of convergence if you're generating
Re: [gmx-users] deltaG from PMF
On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 500 ; kJ mol^-1 nm^-2 pull_nstxout= 1000; every 2 ps pull_nstfout= 1000 ; every 2 ps Based on these settings you're allowing grompp to set the reference distance to whatever it finds in the coordinate file. It seems clear to me that the sampling indicates what I said before - you have an energy minimum somewhere other than where you started with. What that state corresponds to relative to what you think is going on is for you to decide based on the nature of your system. Whatever is occurring at 0.6 nm of COM separation is of particular interest, since the energy minimum is so distinct. So based on this the deltaG will correspond to -5.22 as the initial state was taken at 0.4 nm corresponding to 0 kcal/mol as the moment corresponding to the minimum is the coordinate from SMD where last hydrogen bond was broken. Would you agree? Based on the very little information I have, no. It would appear that the 0.4 nm separation is in fact some metastable state and the true energy minimum is at 0.6 nm of COM separation. What's going on at that location? My mistake. The initial grompp of 1st configuartion (where ligand stakced on keratin surface) corresponds to 0.6 nm where is the minimum. Thus deltaG would be -7.22 kcal/mol. Am I right? Or Shall I take difference between 0 and 5.22 ? -7.22 kcal/mol sounds much more logical to me. If your first configuration is at the energy minimum, that's not something you ignore. The zero point can be set wherever you like with the g_wham flag -zprof0, so it's really rather arbitrary. The WHAM algorithm simply sets the leftmost window (smallest
[gmx-users] deltaG from PMF
Am 21.08.2012 18:22, schrieb gmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull
Re: [gmx-users] deltaG from PMF
Thanks Thomas. Justin, could you please comment on this? Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesier schl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schrieb gmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc
Re: [gmx-users] deltaG from PMF
On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesier schl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schrieb gmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell
