[gmx-users] Re: free energy
Dear Kieu Thu Thanks for your comment about free energy. Unfortunately, I could not send a email to Paissoni Cristina in the Gromacs Forum. Could you give me email address of Paissoni Cristina? Finding a tool for calculation MM/PBSA with Gromacs is very vital for me. Best Regards Kiana -- View this message in context: http://gromacs.5086.x6.nabble.com/free-energy-tp5012246p5012363.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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] Re: free energy calculations of methane in water computed by GMX ver 4.5.7 and ver 4.6.2
On 6/23/13 1:13 AM, Dwey wrote: Hi Justin, Thank you for sharing your experience with me. As suggested, Gromacs ver 4.5.5 is compiled within the same Linux box and I am able to reproduce a similar result ( DG= -9.30 kJmol-1). Gromacs ver 4.5.5 and 4.6.2 both are compiled from source codes, while Gromacs ver 4.5.7 as reported earlier is installed with pre-compiled binary files. Thus, Gromacs ver 4.5.7 is now re-compiled again by myself. Below are shown results. In addition, L-BFGS mdp file works well for all versions after it is modified by adding define = -DFLEXIBLE. Gromacs ver 4.5.5 (compiled from source codes) point 0.050 - 0.100, DG -0.04 +/- 0.00 point 0.100 - 0.150, DG -0.08 +/- 0.00 point 0.150 - 0.200, DG -0.13 +/- 0.01 point 0.200 - 0.250, DG -0.19 +/- 0.00 point 0.250 - 0.300, DG -0.27 +/- 0.00 point 0.300 - 0.350, DG -0.35 +/- 0.00 point 0.350 - 0.400, DG -0.43 +/- 0.01 point 0.400 - 0.450, DG -0.55 +/- 0.01 point 0.450 - 0.500, DG -0.71 +/- 0.01 point 0.500 - 0.550, DG -0.94 +/- 0.00 point 0.550 - 0.600, DG -1.25 +/- 0.00 point 0.600 - 0.650, DG -1.40 +/- 0.01 point 0.650 - 0.700, DG -1.29 +/- 0.01 point 0.700 - 0.750, DG -1.01 +/- 0.00 point 0.750 - 0.800, DG -0.67 +/- 0.00 point 0.800 - 0.850, DG -0.36 +/- 0.00 point 0.850 - 0.900, DG -0.09 +/- 0.00 point 0.900 - 0.950, DG 0.14 +/- 0.00 point 0.950 - 1.000, DG 0.33 +/- 0.00 total 0.050 - 1.000, DG -9.30 +/- 0.03 Gromacs ver 4.5.7 (compiled from source codes) lambda 0.000 - 0.050, DG 0.05 +/- 0.00 lambda 0.050 - 0.100, DG 0.02 +/- 0.00 lambda 0.100 - 0.150, DG -0.04 +/- 0.00 lambda 0.150 - 0.200, DG -0.09 +/- 0.00 lambda 0.200 - 0.250, DG -0.14 +/- 0.01 lambda 0.250 - 0.300, DG -0.21 +/- 0.01 lambda 0.300 - 0.350, DG -0.29 +/- 0.01 lambda 0.350 - 0.400, DG -0.37 +/- 0.00 lambda 0.400 - 0.450, DG -0.48 +/- 0.01 lambda 0.450 - 0.500, DG -0.65 +/- 0.01 lambda 0.500 - 0.550, DG -0.89 +/- 0.01 lambda 0.550 - 0.600, DG -1.19 +/- 0.01 lambda 0.600 - 0.650, DG -1.34 +/- 0.01 lambda 0.650 - 0.700, DG -1.23 +/- 0.00 lambda 0.700 - 0.750, DG -0.95 +/- 0.01 lambda 0.750 - 0.800, DG -0.62 +/- 0.00 lambda 0.800 - 0.850, DG -0.31 +/- 0.00 lambda 0.850 - 0.900, DG -0.03 +/- 0.00 lambda 0.900 - 0.950, DG 0.19 +/- 0.00 lambda 0.950 - 1.000, DG 0.38 +/- 0.00 total 0.000 - 1.000, DG -8.19 +/- 0.03 After comparing the output from ver 4.5.5 with that from ver 4.5.7, I do find quirky information from the g_bar in ver 4.5.7. In ver 4.5.7, for example, it shows that ver 4.5.7 dose not give information of dH/dl (see below) Moreover, I also try g_bar of ver 4.5.5 or 4.6.2 to process the output data (such as md*.xvg generated by ver 4.5.7). The result is unchanged and DG ( -8.19 kJmol-1) remains incorrect. Thanks, Dwey ++ g_bar ver 4.5.7, md0.05.xvg: 0.0 - 5000.0; lambda = 0.050 foreign lambdas: 0.050 (250001 pts) 0.000 (250001 pts) 0.100 (250001 pts) md0.15.xvg: 0.0 - 5000.0; lambda = 0.150 foreign lambdas: 0.150 (250001 pts) 0.100 (250001 pts) 0.200 (250001 pts) md0.1.xvg: 0.0 - 5000.0; lambda = 0.100 foreign lambdas: 0.100 (250001 pts) 0.050 (250001 pts) 0.150 (250001 pts) . . . g_bar ver 4.5.5 or 4.6.2, md0.05.xvg: 0.0 - 5000.0; lambda = 0.05 dH/dl foreign lambdas: dH/dl (250001 pts) delta H to 0 (250001 pts) delta H to 0.1 (250001 pts) md0.15.xvg: 0.0 - 5000.0; lambda = 0.15 dH/dl foreign lambdas: dH/dl (250001 pts) delta H to 0.1 (250001 pts) delta H to 0.2 (250001 pts) md0.1.xvg: 0.0 - 5000.0; lambda = 0.1 dH/dl foreign lambdas: dH/dl (250001 pts) delta H to 0.05 (250001 pts) delta H to 0.15 (250001 pts) . . . Thanks for the detailed information. Please file a bug report on redmine.gromacs.org with all of this information. We had intended for 4.5.7 to be the final version in the 4.5.x series, but if there are serious issues with it, they should probably be corrected. -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
[gmx-users] Re: free energy calculations of methane in water computed by GMX ver 4.5.7 and ver 4.6.2
Hi Justin, Thank you for sharing your experience with me. As suggested, Gromacs ver 4.5.5 is compiled within the same Linux box and I am able to reproduce a similar result ( DG= -9.30 kJmol-1). Gromacs ver 4.5.5 and 4.6.2 both are compiled from source codes, while Gromacs ver 4.5.7 as reported earlier is installed with pre-compiled binary files. Thus, Gromacs ver 4.5.7 is now re-compiled again by myself. Below are shown results. In addition, L-BFGS mdp file works well for all versions after it is modified by adding define = -DFLEXIBLE. Gromacs ver 4.5.5 (compiled from source codes) point 0.050 - 0.100, DG -0.04 +/- 0.00 point 0.100 - 0.150, DG -0.08 +/- 0.00 point 0.150 - 0.200, DG -0.13 +/- 0.01 point 0.200 - 0.250, DG -0.19 +/- 0.00 point 0.250 - 0.300, DG -0.27 +/- 0.00 point 0.300 - 0.350, DG -0.35 +/- 0.00 point 0.350 - 0.400, DG -0.43 +/- 0.01 point 0.400 - 0.450, DG -0.55 +/- 0.01 point 0.450 - 0.500, DG -0.71 +/- 0.01 point 0.500 - 0.550, DG -0.94 +/- 0.00 point 0.550 - 0.600, DG -1.25 +/- 0.00 point 0.600 - 0.650, DG -1.40 +/- 0.01 point 0.650 - 0.700, DG -1.29 +/- 0.01 point 0.700 - 0.750, DG -1.01 +/- 0.00 point 0.750 - 0.800, DG -0.67 +/- 0.00 point 0.800 - 0.850, DG -0.36 +/- 0.00 point 0.850 - 0.900, DG -0.09 +/- 0.00 point 0.900 - 0.950, DG 0.14 +/- 0.00 point 0.950 - 1.000, DG 0.33 +/- 0.00 total 0.050 - 1.000, DG -9.30 +/- 0.03 Gromacs ver 4.5.7 (compiled from source codes) lambda 0.000 - 0.050, DG 0.05 +/- 0.00 lambda 0.050 - 0.100, DG 0.02 +/- 0.00 lambda 0.100 - 0.150, DG -0.04 +/- 0.00 lambda 0.150 - 0.200, DG -0.09 +/- 0.00 lambda 0.200 - 0.250, DG -0.14 +/- 0.01 lambda 0.250 - 0.300, DG -0.21 +/- 0.01 lambda 0.300 - 0.350, DG -0.29 +/- 0.01 lambda 0.350 - 0.400, DG -0.37 +/- 0.00 lambda 0.400 - 0.450, DG -0.48 +/- 0.01 lambda 0.450 - 0.500, DG -0.65 +/- 0.01 lambda 0.500 - 0.550, DG -0.89 +/- 0.01 lambda 0.550 - 0.600, DG -1.19 +/- 0.01 lambda 0.600 - 0.650, DG -1.34 +/- 0.01 lambda 0.650 - 0.700, DG -1.23 +/- 0.00 lambda 0.700 - 0.750, DG -0.95 +/- 0.01 lambda 0.750 - 0.800, DG -0.62 +/- 0.00 lambda 0.800 - 0.850, DG -0.31 +/- 0.00 lambda 0.850 - 0.900, DG -0.03 +/- 0.00 lambda 0.900 - 0.950, DG 0.19 +/- 0.00 lambda 0.950 - 1.000, DG 0.38 +/- 0.00 total 0.000 - 1.000, DG -8.19 +/- 0.03 After comparing the output from ver 4.5.5 with that from ver 4.5.7, I do find quirky information from the g_bar in ver 4.5.7. In ver 4.5.7, for example, it shows that ver 4.5.7 dose not give information of dH/dl (see below) Moreover, I also try g_bar of ver 4.5.5 or 4.6.2 to process the output data (such as md*.xvg generated by ver 4.5.7). The result is unchanged and DG ( -8.19 kJmol-1) remains incorrect. Thanks, Dwey ++ g_bar ver 4.5.7, md0.05.xvg: 0.0 - 5000.0; lambda = 0.050 foreign lambdas: 0.050 (250001 pts) 0.000 (250001 pts) 0.100 (250001 pts) md0.15.xvg: 0.0 - 5000.0; lambda = 0.150 foreign lambdas: 0.150 (250001 pts) 0.100 (250001 pts) 0.200 (250001 pts) md0.1.xvg: 0.0 - 5000.0; lambda = 0.100 foreign lambdas: 0.100 (250001 pts) 0.050 (250001 pts) 0.150 (250001 pts) . . . g_bar ver 4.5.5 or 4.6.2, md0.05.xvg: 0.0 - 5000.0; lambda = 0.05 dH/dl foreign lambdas: dH/dl (250001 pts) delta H to 0 (250001 pts) delta H to 0.1 (250001 pts) md0.15.xvg: 0.0 - 5000.0; lambda = 0.15 dH/dl foreign lambdas: dH/dl (250001 pts) delta H to 0.1 (250001 pts) delta H to 0.2 (250001 pts) md0.1.xvg: 0.0 - 5000.0; lambda = 0.1 dH/dl foreign lambdas: dH/dl (250001 pts) delta H to 0.05 (250001 pts) delta H to 0.15 (250001 pts) . . . On 6/21/13 11:07 AM, Dwey wrote: Hi gmx-users, I almost reproduced free energy calculations of methane in water on Justin's website. First of all, I am able to follow the workflow of computing solvation free energy for several times with Gromacs version 4.5.7 and version 4.6.2 installed in two identical Linux boxes. However. the output results of GMX ver 4.5.7 and ver 4.6.2 show different values of dG ## GMX Ver. 4.5.7: lambda 0.000 - 0.050, DG 0.05 +/- 0.00 lambda 0.050 - 0.100, DG 0.01 +/- 0.00 lambda 0.100 - 0.150, DG -0.03 +/- 0.01 lambda 0.150 - 0.200, DG -0.08 +/- 0.00 lambda 0.200 - 0.250, DG -0.15 +/- 0.00 lambda 0.250 - 0.300, DG -0.21 +/- 0.01 lambda 0.300 - 0.350, DG -0.28 +/- 0.00 lambda 0.350 - 0.400, DG -0.38 +/- 0.00 lambda 0.400 - 0.450,
[gmx-users] Re: Free Energy Calculations in Gromacs
Thanks Justin and Professor Shirts. The gmxdump on the resulting tpr did indeed solve a lot of my problems. Thanks again. -- View this message in context: http://gromacs.5086.x6.nabble.com/Free-Energy-Calculations-in-Gromacs-tp5007495p5009115.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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] Re: Free Energy Calculations in Gromacs
Hi, I have been trying to perform the simulations using the amber forcefield, in which the [ pairtypes ] directive is not defined explicitly in the ffnonbonded.itp file, but rather the 1-4 interactions are generated as per the defaults section in the forcefield.itp file. I was wondering what happens to these 1-4 interactions at lambda=1 state? Suppose 1 corresponds to CA and 4 corresponds to NA at state A (lambda = 0) and the CA_per and NA_per are the corresponding atomtypes at state B (lambda = 1). It is defined in the topology file that ... ... [ pairs ] 1 4 1 So does it mean that at state B (lambda = 1), the 1-4 nonbonded interactions will be calculated between CA_per and NA_per? Although the [ pairs_nb ] section is described in the manual, I would appreciate if someone elaborates a little more. -- View this message in context: http://gromacs.5086.x6.nabble.com/Free-Energy-Calculations-in-Gromacs-tp5007495p5008996.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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] Re: Free Energy Calculations in Gromacs
On 6/10/13 2:50 PM, JW Gibbs wrote: Hi, I have been trying to perform the simulations using the amber forcefield, in which the [ pairtypes ] directive is not defined explicitly in the ffnonbonded.itp file, but rather the 1-4 interactions are generated as per the defaults section in the forcefield.itp file. I was wondering what happens to these 1-4 interactions at lambda=1 state? Suppose 1 corresponds to CA and 4 corresponds to NA at state A (lambda = 0) and the CA_per and NA_per are the corresponding atomtypes at state B (lambda = 1). It is defined in the topology file that ... ... [ pairs ] 1 4 1 So does it mean that at state B (lambda = 1), the 1-4 nonbonded interactions will be calculated between CA_per and NA_per? Although the [ pairs_nb ] section is described in the manual, I would appreciate if someone elaborates a little more. The actual answer depends on exactly how you're treating the system. Are you using [pairs_nb] in your topology or just [pairs]? The outcome will be different depending on which you're using. Also note that, as the manual says, you don't really need to mess with [pairs_nb] at all; you can achieve unscaled internal interactions using couple-intramol = no in the .mdp file. Without seeing the .mdp file, it's even more difficult to know what you're doing and what you should expect. -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
Re: [gmx-users] Re: Free Energy Calculations in Gromacs
An important final point is that you can always see EXACTLY what grompp is putting into the B state by running gmxdump on the resulting tpr. It's a LOT of information, but all in text all the interactions are listed explicitly there. On Mon, Jun 10, 2013 at 6:20 PM, Justin Lemkul jalem...@vt.edu wrote: On 6/10/13 2:50 PM, JW Gibbs wrote: Hi, I have been trying to perform the simulations using the amber forcefield, in which the [ pairtypes ] directive is not defined explicitly in the ffnonbonded.itp file, but rather the 1-4 interactions are generated as per the defaults section in the forcefield.itp file. I was wondering what happens to these 1-4 interactions at lambda=1 state? Suppose 1 corresponds to CA and 4 corresponds to NA at state A (lambda = 0) and the CA_per and NA_per are the corresponding atomtypes at state B (lambda = 1). It is defined in the topology file that ... ... [ pairs ] 1 4 1 So does it mean that at state B (lambda = 1), the 1-4 nonbonded interactions will be calculated between CA_per and NA_per? Although the [ pairs_nb ] section is described in the manual, I would appreciate if someone elaborates a little more. The actual answer depends on exactly how you're treating the system. Are you using [pairs_nb] in your topology or just [pairs]? The outcome will be different depending on which you're using. Also note that, as the manual says, you don't really need to mess with [pairs_nb] at all; you can achieve unscaled internal interactions using couple-intramol = no in the .mdp file. Without seeing the .mdp file, it's even more difficult to know what you're doing and what you should expect. -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] Re: Free energy of solvation about peptide
The BAR method calculates a ratio of weighted average of the Hamiltonian difference of state B given state A and vice versa. So, for different peptides, can this method apply to compare the stability of peptides? Thank you very much! -- View this message in context: http://gromacs.5086.x6.nabble.com/Free-energy-of-solvation-about-peptide-tp5008353p5008437.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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] Re: free energy of system
Dear Justin, Thanks for your response and the perfect tutorial. I found the details I needed in it. Best regards, Sarah -- 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] Re: Free Energy tutorial - choosing number of solvent molecules
Hello Justin, I am running 2,000,000 time steps for the actual MD run (4,000 ps). Depending on how many solvent molecules I start with, I get slightly different results. Do you think its ok to run several different tests and take the average? Or perhaps take the end results of a shorter MD run and use those as the starting coordinates for a new run? Thanks, Fabian Fabian Casteblanco wrote: Hello all, I'm running the same process from the free energy tutorial by Justin Lemkul...http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html. How did the number of solvent particles get chosen (in the tutorial, 210 molecules were chosen)? I seem to be getting slightly different If memory serves, I reproduced what was in the original paper, but do check. results (ranging from as small as 200 J/mol to about 5 kJ/mol depending on how many molecules I choose (ranging for example from 210 ethanol molecules to about 610 ethanol molecules for the largest energy difference change of about 5 kJ). I keep running tests to see if there is some sort of minimum atom number to get steady consistent numbers but I can't seem to find it. When I plot for example the bar.xvg barint.xvg for both sets to see where the lines don't match up, its usually one or two points that differ slightly which cause the free energies in the end to be slightly different.I seem to be noticing too that the more atoms I use, the free energy gets a little bit lower. Does anybody have any experience with this? How long are your simulations? I have experienced the case (using a water-octanol solvent mixture) where the initial configuration made a big difference in the result, so longer simulations and multiple configurations for the solvent were necessary to get reliable averages. -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin On Thu, Jan 26, 2012 at 11:21 AM, Fabian Casteblanco fabian.castebla...@gmail.com wrote: Hello all, I'm running the same process from the free energy tutorial by Justin Lemkul...http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html. How did the number of solvent particles get chosen (in the tutorial, 210 molecules were chosen)? I seem to be getting slightly different results (ranging from as small as 200 J/mol to about 5 kJ/mol depending on how many molecules I choose (ranging for example from 210 ethanol molecules to about 610 ethanol molecules for the largest energy difference change of about 5 kJ). I keep running tests to see if there is some sort of minimum atom number to get steady consistent numbers but I can't seem to find it. When I plot for example the bar.xvg barint.xvg for both sets to see where the lines don't match up, its usually one or two points that differ slightly which cause the free energies in the end to be slightly different. I seem to be noticing too that the more atoms I use, the free energy gets a little bit lower. Does anybody have any experience with this? Thanks. -- Best regards, Fabian F. Casteblanco Rutgers University -- C: +908 917 0723 E: fabian.castebla...@gmail.com -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Re: Free Energy tutorial - choosing number of solvent molecules
Fabian Casteblanco wrote: Hello Justin, I am running 2,000,000 time steps for the actual MD run (4,000 ps). Depending on how many solvent molecules I start with, I get slightly different results. Do you think its ok to run several different tests and take the average? Or perhaps take the end results of a shorter MD run and use those as the starting coordinates for a new run? Either approach sounds reasonable to me. Keep in mind magnitude - some of your differences were 5 kJ/mol, which may or may not be statistically relevant, depending on the value of DG and the resulting error estimates. -Justin Thanks, Fabian Fabian Casteblanco wrote: Hello all, I'm running the same process from the free energy tutorial by Justin Lemkul...http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html. How did the number of solvent particles get chosen (in the tutorial, 210 molecules were chosen)? I seem to be getting slightly different If memory serves, I reproduced what was in the original paper, but do check. results (ranging from as small as 200 J/mol to about 5 kJ/mol depending on how many molecules I choose (ranging for example from 210 ethanol molecules to about 610 ethanol molecules for the largest energy difference change of about 5 kJ). I keep running tests to see if there is some sort of minimum atom number to get steady consistent numbers but I can't seem to find it. When I plot for example the bar.xvg barint.xvg for both sets to see where the lines don't match up, its usually one or two points that differ slightly which cause the free energies in the end to be slightly different.I seem to be noticing too that the more atoms I use, the free energy gets a little bit lower. Does anybody have any experience with this? How long are your simulations? I have experienced the case (using a water-octanol solvent mixture) where the initial configuration made a big difference in the result, so longer simulations and multiple configurations for the solvent were necessary to get reliable averages. -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee 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] Re: Free energy calculation question
Hello Justin, I'm looking at the dG vs Lambda plot that is an output from G_bar but on the Shirts et al paper that you included, the part where it talks about linearity, it is referring to dH/dLambda for electrostatic decoupling. So if I take the line formed by dGtotal vs Lambda from g_bar output and take the derivative excluding the beginning and end parts (lambda approaches 0 and 1) then it does seem to take more of a somewhat linear shape. If you can see attached, I think perhaps thats how I was misreading it? I hope that can be the reason, since I've done it with 4 different cases (2 different amount of solvent molecules and using MD and SD integrator) and they all seem to be giving similar results. Thanks for your help. -Fabian On Wed, Aug 31, 2011 at 11:24 AM, Fabian Casteblanco fabian.castebla...@gmail.com wrote: Hello Justin, I'm calculating the free energy of a drug in an alcohol solvent. I have a question referring to your free energy tutorial. You mentioned that decoupling of electrostatic interactions is linear and decoupling of vdW can vary. Is this true for your case of methanol in water or for all cases? When I ran my system of drug in ethanol solvent, I got a non linear dG for both electrostatic and vdW. Also, is there no need to find dG of cav ( the free energy required to form the solute cavity within the solvent) ? I have attached some graphs. Thanks for your help. Your tutorial was extremely useful. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Re: Free energy calculation question
Fabian Casteblanco wrote: Hello Justin, I'm looking at the dG vs Lambda plot that is an output from G_bar but on the Shirts et al paper that you included, the part where it talks about linearity, it is referring to dH/dLambda for electrostatic decoupling. So if I take the line formed by dGtotal vs Lambda from g_bar output and take the derivative excluding the beginning and end parts (lambda approaches 0 and 1) then it does seem to take more of a somewhat linear shape. If you can see attached, I think perhaps thats You haven't attached anything. how I was misreading it? I hope that can be the reason, since I've done it with 4 different cases (2 different amount of solvent molecules and using MD and SD integrator) and they all seem to be giving similar results. That could be it. The BAR output and the plots from TI are different. Also be sure that you're not using soft-core for the Coulombic (de)coupling, as was suggested previously. -Justin Thanks for your help. -Fabian On Wed, Aug 31, 2011 at 11:24 AM, Fabian Casteblanco fabian.castebla...@gmail.com wrote: Hello Justin, I'm calculating the free energy of a drug in an alcohol solvent. I have a question referring to your free energy tutorial. You mentioned that decoupling of electrostatic interactions is linear and decoupling of vdW can vary. Is this true for your case of methanol in water or for all cases? When I ran my system of drug in ethanol solvent, I got a non linear dG for both electrostatic and vdW. Also, is there no need to find dG of cav ( the free energy required to form the solute cavity within the solvent) ? I have attached some graphs. Thanks for your help. Your tutorial was extremely useful. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee 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] Re: Free energy calculation question
Hello Justin, You mean that only for vdW decoupling, you would need to use soft-core potentials? I had soft -core potentials on for decoupling the electrostatic interactions (see below). What would I use in its place? Thanks again for your help! ;Production MD ; include =-I/mphase/users2/fabian/CGenff title =CGenFF Lov/Eth Solution MD run ;PARAMETERS - describing what to do, when to stop and what to save ; ;Run parameters integrator =sd ;leap-frog integrator ld_seed =-1 nsteps =50 ;2*50=1000ps, 1 ns dt =0.002 ;2 fs nstcomm = 100 ;*** - frequency for center of mass motion removal ;Output control nstxout =1000 ;save coordinates every 2ps nstvout =1000 ;save velocities every 2 ps nstenergy =1000 ;save energies every 2 ps nstlog =1000 ;update log file every 2 ps nstxtcout =1000 ;xtc compressed trajectory output every 2 ps xtc-precision =1000 ;*** - precision to write to xtc trajectory ;Bond Parameters continuation=yes;Restarting after NPT constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds) constr;ained lincs_iter =1 ;accuracy of LINCS lincs_order =12 ;also related to accuracy ;Neighborhood searching ns_type =grid ;search neighboring grid cells nstlist =5 ;#steps. 5*0.002 ps = 5* 2 fs = 10 fs - Frequency to update the neighbor list (and the long-range forces, when ;using twin-range cut-off’s). When this is 0, the neighbor list is made only once. rlist =1.1;short-range neighborlist cutoff (in nm) rcoulomb=1.1;short-range electrostatic cutoff (in nm) pbc =xyz; 3-D PBC ;Electrostatics coulombtype =PME;Particle Mesh Ewald for long-range electrostat;ics pme_order =4 ;cubic interpolation fourierspacing =0.16 ;grid spacing for FFT ; van der Waals vdwtype =Shift ;Van der Waals for CHARMM rvdw_switch =0.8 rvdw=1.0;Short-range Van der Waals cut-off ;Dispersion correction DispCorr=EnerPres ;account for cut-off vdW scheme ;Temperature coupling is on tcoupl =V-rescale ;modified Berendsen thermostat tc-grps =SYSTEM ;two coupling groups - more accurate tau_t =0.1;time constant, in ps ref_t =298;reference temperature, on for each group, in K ;Pressure coupling is on pcoupl=Parrinello-Rahman;Pressure coupling on in NPT pcoupltype =isotropic ;uniform scaling of box vect;ors tau_p =2.0 ;time constant, in ps ref_p =1.0 ;reference pressure, in bar compressibility =4.5e-5;isothermal compr of H2O, ba;r^(-1) ; Free energy control stuff free_energy = yes ;*** - Indicates we are doing a free energy calculation, and that interpolation between the A and B states of the ;chosen molecule (defined below) will occur. init_lambda = 0.0 ;*** - Value of λ delta_lambda= 0 ;*** - The value of λ can be incremented by some amount per timestep (i.e., δλ/δt) in a technique called slow ;growth. This method can have significant errors associated with it, and thus we will make no time-dependent ;changes to our λ values. foreign_lambda = 0.05 ;*** - Additional values of λ for which ΔH will be written to dhdl.xvg (with frequency nstdhdl). The ;configurations generated in the trajectory at λ = init_lambda will have ΔH calculated for these same ;configurations at all values of λ = foreign_lambda. sc-alpha= 0.5 ;*** - the soft-core parameter, a value of 0 results in linear interpolation of the LJ and Coulomb ;interactions sc-power= 1.0 ;*** - the power for lambda in the soft-core function, only the values 1 and 2 are supported sc-sigma
Re: [gmx-users] Re: Free energy calculation question
Fabian Casteblanco wrote: Hello Justin, You mean that only for vdW decoupling, you would need to use soft-core potentials? I had soft -core potentials on for decoupling the electrostatic interactions (see below). What would I use in its place? Thanks again for your help! You only need soft-core for LJ transformations. Set sc_alpha to 0 (as you have noted in the comment) for normal (linear) interpolation. -Justin ;Production MD ; include =-I/mphase/users2/fabian/CGenff title =CGenFF Lov/Eth Solution MD run ;PARAMETERS - describing what to do, when to stop and what to save ; ;Run parameters integrator =sd ;leap-frog integrator ld_seed =-1 nsteps =50 ;2*50=1000ps, 1 ns dt =0.002 ;2 fs nstcomm = 100 ;*** - frequency for center of mass motion removal ;Output control nstxout =1000 ;save coordinates every 2ps nstvout =1000 ;save velocities every 2 ps nstenergy =1000 ;save energies every 2 ps nstlog =1000 ;update log file every 2 ps nstxtcout =1000 ;xtc compressed trajectory output every 2 ps xtc-precision =1000 ;*** - precision to write to xtc trajectory ;Bond Parameters continuation=yes;Restarting after NPT constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds) constr;ained lincs_iter =1 ;accuracy of LINCS lincs_order =12 ;also related to accuracy ;Neighborhood searching ns_type =grid ;search neighboring grid cells nstlist =5 ;#steps. 5*0.002 ps = 5* 2 fs = 10 fs - Frequency to update the neighbor list (and the long-range forces, when ;using twin-range cut-off’s). When this is 0, the neighbor list is made only once. rlist =1.1;short-range neighborlist cutoff (in nm) rcoulomb=1.1;short-range electrostatic cutoff (in nm) pbc =xyz; 3-D PBC ;Electrostatics coulombtype =PME;Particle Mesh Ewald for long-range electrostat;ics pme_order =4 ;cubic interpolation fourierspacing =0.16 ;grid spacing for FFT ; van der Waals vdwtype =Shift ;Van der Waals for CHARMM rvdw_switch =0.8 rvdw=1.0;Short-range Van der Waals cut-off ;Dispersion correction DispCorr=EnerPres ;account for cut-off vdW scheme ;Temperature coupling is on tcoupl =V-rescale ;modified Berendsen thermostat tc-grps =SYSTEM ;two coupling groups - more accurate tau_t =0.1;time constant, in ps ref_t =298;reference temperature, on for each group, in K ;Pressure coupling is on pcoupl=Parrinello-Rahman;Pressure coupling on in NPT pcoupltype =isotropic ;uniform scaling of box vect;ors tau_p =2.0 ;time constant, in ps ref_p =1.0 ;reference pressure, in bar compressibility =4.5e-5;isothermal compr of H2O, ba;r^(-1) ; Free energy control stuff free_energy = yes ;*** - Indicates we are doing a free energy calculation, and that interpolation between the A and B states of the ;chosen molecule (defined below) will occur. init_lambda = 0.0 ;*** - Value of λ delta_lambda= 0 ;*** - The value of λ can be incremented by some amount per timestep (i.e., δλ/δt) in a technique called slow ;growth. This method can have significant errors associated with it, and thus we will make no time-dependent ;changes to our λ values. foreign_lambda = 0.05 ;*** - Additional values of λ for which ΔH will be written to dhdl.xvg (with frequency nstdhdl). The ;configurations generated in the trajectory at λ = init_lambda will have ΔH calculated for these same ;configurations at all values of λ = foreign_lambda. sc-alpha= 0.5 ;*** - the soft-core parameter, a value of 0 results in linear interpolation of the LJ and Coulomb
[gmx-users] Re: Free energy calculation question
Thanks Justin! On Wed, Aug 31, 2011 at 11:24 AM, Fabian Casteblanco fabian.castebla...@gmail.com wrote: Hello Justin, I'm calculating the free energy of a drug in an alcohol solvent. I have a question referring to your free energy tutorial. You mentioned that decoupling of electrostatic interactions is linear and decoupling of vdW can vary. Is this true for your case of methanol in water or for all cases? When I ran my system of drug in ethanol solvent, I got a non linear dG for both electrostatic and vdW. Also, is there no need to find dG of cav ( the free energy required to form the solute cavity within the solvent) ? I have attached some graphs. Thanks for your help. Your tutorial was extremely useful. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Re: Free energy calculation
Thanks Justin. Best regards, Fabian On Mon, Aug 8, 2011 at 5:49 PM, Fabian Casteblanco fabian.castebla...@gmail.com wrote: Hello all, I am setting up a free energy calculation (drug from full coulomb+vdW in solution -- drug with only vdW in solution -- dummy drug in solution). After reading most of the papers, I understand that you need significant overlap from the energies for each intermediate point to overlap so its best to have many intermediate points from Lambda=0 to Lambda=1. The drug molecule I have is a bit complex but I wasn't sure if using too small of an intermediate could have a bad effect on the free_energy calculation. I know Justin Lemkul said in its tutorial that Lambda=+0.05 should be good for most but I decided to go with Lambda=+0.02. Could this have any negative effect other than taking a longer time? Also, how does one come up with the best soft-core potential parameters to use? Is it ok to use the one from the Methane in water tutorial? Thanks for everyone's expertise. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Re: free energy calculations
Message: 1 Date: Tue, 05 Jul 2011 19:14:04 -0400 From: Justin A. Lemkul jalem...@vt.edu Subject: Re: [gmx-users] free energy calculations To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: 4e139abc.1030...@vt.edu Content-Type: text/plain; charset=ISO-8859-1; format=flowed Ragothaman Yennamalli wrote: Hi, I want to run free energy calculations on a particular protein-ligand complex. I do not have much knowledge on this so I have some questions, hopefully someone might give me clear answers. I am following the tutorial given in this link: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/08_advanced.html My questions are as follows: 1. I am assuming that I do not need the protein without the ligand form for free energy calculations. Am I right? Free energy calculations can be done on a number of systems. The tutorial is just one simple example. Thanks for the mail. So, I technically do not need both the un-complexed and complexed form of the protein with the ligand. 2. To use g_bar, I need to run the protein-ligand complex with a) first with lambda ranging from 0.5 to 1 for vdw coupling and b) followed by lambda ranging from 0.5 to 1 for coulombic coupling Adding these would give the deltaG I am looking for. Am I right? No. The range for lambda should be 0 to 1. Using 0.5 to 1 gives only part of the free energy for the transformation, and not necessarily half. Apologies, I should have mentioned as from 0 to 1. So, the first step is to do vdw coupling followed by coulombic coupling. Right? 3. Are there any papers that use g_bar function in gromacs to calculate free energy? The g_bar tool is relatively new to Gromacs, so maybe not. The BAR method itself has been around for decades though, so its applications have been demonstrated in the literature many times. Thanks for pointing this. Raghu -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee 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
Re: [gmx-users] Re: free energy calculations
Ragothaman Yennamalli wrote: Message: 1 Date: Tue, 05 Jul 2011 19:14:04 -0400 From: Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu Subject: Re: [gmx-users] free energy calculations To: Discussion list for GROMACS users gmx-users@gromacs.org mailto:gmx-users@gromacs.org Message-ID: 4e139abc.1030...@vt.edu mailto:4e139abc.1030...@vt.edu Content-Type: text/plain; charset=ISO-8859-1; format=flowed Ragothaman Yennamalli wrote: Hi, I want to run free energy calculations on a particular protein-ligand complex. I do not have much knowledge on this so I have some questions, hopefully someone might give me clear answers. I am following the tutorial given in this link: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/08_advanced.html My questions are as follows: 1. I am assuming that I do not need the protein without the ligand form for free energy calculations. Am I right? Free energy calculations can be done on a number of systems. The tutorial is just one simple example. Thanks for the mail. So, I technically do not need both the un-complexed and complexed form of the protein with the ligand. This is discussed in the tutorial: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/08_advanced.html 2. To use g_bar, I need to run the protein-ligand complex with a) first with lambda ranging from 0.5 to 1 for vdw coupling and b) followed by lambda ranging from 0.5 to 1 for coulombic coupling Adding these would give the deltaG I am looking for. Am I right? No. The range for lambda should be 0 to 1. Using 0.5 to 1 gives only part of the free energy for the transformation, and not necessarily half. Apologies, I should have mentioned as from 0 to 1. So, the first step is to do vdw coupling followed by coulombic coupling. Right? No, the opposite. You do not want to remove vdW terms from charged species. This topic is addressed in the tutorial. http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/02_topology.html -Justin 3. Are there any papers that use g_bar function in gromacs to calculate free energy? The g_bar tool is relatively new to Gromacs, so maybe not. The BAR method itself has been around for decades though, so its applications have been demonstrated in the literature many times. Thanks for pointing this. Raghu -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu http://vt.edu | (540) 231-9080 tel:%28540%29%20231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee 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] Re: free energy
Moeed, Sorry for the late response, and yes I did need to include more information in my answer. 1) If your polymers are physically close together in reality you might want to have multiple polymers in a single FEP job. I run single ligands in proteins, so this has never been an issue. You may want to run both and see what the differences in energy are. If there are none you may be able to do longer md/sd runs with a smaller one polymer system than you could with multiple polymers in the same time. In other words I have no idea what your system will need. Sorry. 2) How many windows and how long are the md/sd runs for each window? I personally found that 21 windows works well for my systems. This however will vary by the dV/dlambda graph. If and when your dV/dlambda is changing more (ddV/dlambda) you will need more windows to get a more accurate energy. Most systems change the most at the beginning and end of lambda values. I run a lambda every 0.05, but my systems are huge! (100,000 atoms) 3) Bennett acceptance ratio is an extra step in setup but makes getting your data much easier to interpret. Once you are done with all your windows you run g_bar and it will print out on your screen the energy from each window and the total energy in kJ/mol with a standard deviation. This can be found in the manual and also in the gmx_user archives. This entailes making a .mdp file for every step (but you probably do this already) and setting the foreign lambdas correctly. ie for lambda values: 0.00, 0.25, 0.50, 0.75 and 1.00 Your foreign lambda values for the 0.00 mdp file would be: 0.25 Your foreign lambda values for the 0.25 mdp file would be: 0.00 0.50 Your foreign lambda values for the 0.50 mdp file would be: 0.25 0.75 Your foreign lambda values for the 0.75 mdp file would be: 0.50 1.00 Your foreign lambda values for the 0.50 mdp file would be: 0.75 This of course is a small number of lambda values you will want to use more. Please email me with questions if you need more. Thank you, TJ Mustard Email: musta...@onid.orst.edu PS more below On February 14, 2011 at 8:54 PM Moeed lecie...@googlemail.com wrote: Hello, Sorry but it seems you wanted to answer my question on gmx list but I dont see anything!... -- Moeed Shahamat Graduate Student (Materials Modeling Research Group) McGill University- Department of Chemical Engineering Montreal, Quebec H3A 2B2, Canada Web:http://mmrg.chemeng.mcgill.ca/pages/current-group-members/moeed-shahamat.php Web:http://mmrg.chemeng.mcgill.ca/ TJ Mustard Email: musta...@onid.orst.edu Dear experts, I am going to do solvation FE of polymer (polyethylene) in a hydrocarbon solvent. I have prepared a system consisting of 4 polymer chains and 480 hexane molecules with the actual density of polymer solution (~ 0.5 g/cm3). 1- For such a study I dont know how many polymers I need to have in my system. If FE can be done with only one chain, am I making system bigger in vain? Does this matter affect the accuracy of results? 2- I have switched off electrostatics so I am using free_energy = yes init_lambda = 0 delta_lambda = 0 sc_alpha = 0.5 sc-power = 1 sc_sigma = 0.3 couple-lambda0 = vdw couple-lambda1 = none couple-intramol = no In David Mobleys turorial the last three lines are not included. I wanted to know if I am to run say 10 simulations for different lambda, what purpose does the last three lines serve in 4.0.7 ? I got very close values in that tutorial without these settings. ( I know what these lines mean, just curious how these three lines affect the results in 4 X +). This old tutorial is for an old gromacs which needs you to set the state a and state b. In the newer gromacs releases you set the state a, and gromacs will go to state b without you setting each atom. This considerably saves time making .itp files. All you need is a methane .itp file (got mine from antechamber) and then you set couple-lambda0 and couple-lambda1 to what you are looking for. I personally go from couple-lambda0 = vdw-q to couple-lambda1 = none. Please let me know your comments/point of view about the system and setting I am using. Thanks Moeed TJ Mustard Email: musta...@onid.orst.edu
Re: [gmx-users] Re: free energy calculation , grompp crash
hello Now it is almost clear what happened. When couple-intramol is no (default), all pairwise vdm and charge interaction becomes bonded interaction. All intra-molecular non-bonded interactions for moleculetype couple-moltype are replaced by exclusions and explicit pair interactions. In this manner the decoupled state of the molecule corresponds to the proper vacuum state without periodicity effects. But there is still one thing I don't fully understand. Mdrun will complain Warning: 1-4 interaction between 1 and 114 at distance 2.035 which is larger than the 1-4 table size 2.000 nm. I guess is that every pair wiii have 1-4 interaction because every pair (except exclusion) becomes bonded interaction. These fake 1-4 are not calculated as an 1-4 interaction so that I can neglect this warning, right? thanks. dawei -- 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] Re: free energy calculation , grompp crash
hi, I did more test and found that it depended on size of the protein. Grompp will die when number of atoms of the protein is larger than about 200. Is it possible the source code limit the size of the protein that can be decoupled? thanks. dawei On Mon, Feb 7, 2011 at 10:34 AM, Da-Wei Li lida...@gmail.com wrote: Dear users I tried free energy calculation but grompp couldn't go through. It stops after *** Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein' turning H bonds into constraints... Excluding 2 bonded neighbours molecule type 'SOL' turning H bonds into constraints... Excluding 1 bonded neighbours molecule type 'CL' turning H bonds into constraints... Coupling 1 copies of molecule type 'Protein' *** The CPU usage is 100%. I just add following into the mdp file: *** free_energy = yes init_lambda = 0.0 delta_lambda = 0 sc_alpha =0.5 sc-power =1.0 sc-sigma = 0.3 couple-moltype = Protein couple-lambda0 = vdw-q couple-lambda1 = none *** Does anyone have some idea about this problem? thanks. Another question is whether I can switch off two molecules (such as protein+ligand) in free energy calculation? I searched this list and got that 4.0.7 did support this. how about 4.5.4? dawei -- 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] Re: free energy calculation , grompp crash
Well. It actually isn't dead but becomes very slow for large proteins. dawei On Mon, Feb 7, 2011 at 11:44 AM, Da-Wei Li lida...@gmail.com wrote: hi, I did more test and found that it depended on size of the protein. Grompp will die when number of atoms of the protein is larger than about 200. Is it possible the source code limit the size of the protein that can be decoupled? thanks. dawei On Mon, Feb 7, 2011 at 10:34 AM, Da-Wei Li lida...@gmail.com wrote: Dear users I tried free energy calculation but grompp couldn't go through. It stops after *** Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein' turning H bonds into constraints... Excluding 2 bonded neighbours molecule type 'SOL' turning H bonds into constraints... Excluding 1 bonded neighbours molecule type 'CL' turning H bonds into constraints... Coupling 1 copies of molecule type 'Protein' *** The CPU usage is 100%. I just add following into the mdp file: *** free_energy = yes init_lambda = 0.0 delta_lambda = 0 sc_alpha =0.5 sc-power =1.0 sc-sigma = 0.3 couple-moltype = Protein couple-lambda0 = vdw-q couple-lambda1 = none *** Does anyone have some idea about this problem? thanks. Another question is whether I can switch off two molecules (such as protein+ligand) in free energy calculation? I searched this list and got that 4.0.7 did support this. how about 4.5.4? dawei -- 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] Re: free energy calculation , grompp crash
Dawei, I have no problems with proteins in the thousands of atoms. Can you post your command line and mdp files? Thank you, TJ Mustard On February 7, 2011 at 9:31 AM Da-Wei Li lida...@gmail.com wrote: Well. It actually isnt dead but becomes very slow for large proteins. dawei On Mon, Feb 7, 2011 at 11:44 AM, Da-Wei Li lida...@gmail.com wrote: hi, I did more test and found that it depended on size of the protein. Grompp will die when number of atoms of the protein is larger than about 200. Is it possible the source code limit the size of the protein that can be decoupled? thanks. dawei On Mon, Feb 7, 2011 at 10:34 AM, Da-Wei Li lida...@gmail.com wrote: Dear users I tried free energy calculation but grompp couldnt go through. It stops after *** Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type Protein turning H bonds into constraints... Excluding 2 bonded neighbours molecule type SOL turning H bonds into constraints... Excluding 1 bonded neighbours molecule type CL turning H bonds into constraints... Coupling 1 copies of molecule type Protein *** The CPU usage is 100%. I just add following into the mdp file: *** free_energy = yes init_lambda = 0.0 delta_lambda = 0 sc_alpha =0.5 sc-power =1.0 sc-sigma = 0.3 couple-moltype = Protein couple-lambda0 = vdw-q couple-lambda1 = none *** Does anyone have some idea about this problem? thanks. Another question is whether I can switch off two molecules (such as protein+ligand) in free energy calculation? I searched this list and got that 4.0.7 did support this. how about 4.5.4? dawei TJ Mustard Email: musta...@onid.orst.edu -- 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
Re: [gmx-users] Re: free energy calculation , grompp crash
Da-Wei Li wrote: hello Here they are the command line and mdp file. I use Gromacs 4.5.3. This is a test case only and the protein is 1UBQ. Grompp wills top for about 10 minutes then go through. The efficiency of this kind of process will depend on the amount of available memory on the system. You're asking grompp to decouple a huge amount of degrees of freedom, which will require a lot of memory to do. It sounds like it's working, in any case, so there's no real problem. Whether or not simultaneously decoupling the LJ and Coulombic interactions of a whole protein will generate a stable trajectory or sensible result is another matter. -Justin ***output of grompp* Back Off! I just backed up mdout.mdp to ./#mdout.mdp.6# Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein' turning H bonds into constraints... Excluding 2 bonded neighbours molecule type 'SOL' turning H bonds into constraints... Coupling 1 copies of molecule type 'Protein' Setting gen_seed to 8552 Velocities were taken from a Maxwell distribution at 300 K Command line and mdp file: ** grompp -f pr1.mdp -c after_em.gro -t em.trr -p topol.top -o pr1 ** define = -DPOSRES ; position restrain the protein ; Run parameters integrator = sd ; leap-frog integrator nsteps = 5000 ; 2 * 5 = 100 ps dt = 0.002 ; 2 fs ; Output control nstxout = 1000 ; save coordinates every 2 ps nstvout = 5000 ; save velocities every 100ps nstenergy = 1000 ; save energies every 2 ps nstlog = 1000 ; update log file every 2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = hbonds ; H bonds constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cels nstlist = 10 ; 20 fs rlist = 0.8 ; short-range neighborlist cutoff (in nm) rcoulomb = 0.8 ; short-range electrostatic cutoff (in nm) rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein Non-Protein ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no ; no pressure coupling in NVT ; 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 = 300 ; temperature for Maxwell distribution gen_seed = -1 ; generate a random seed ;free energy stuff free_energy = yes init_lambda = 0.0 delta_lambda = 0 sc_alpha =0.5 sc-power =1.0 sc-sigma = 0.3 couple-moltype = Protein couple-lambda0 = vdw-q couple-lambda1 = none *** thanks. dawei On Mon, Feb 7, 2011 at 1:05 PM, TJ Mustard musta...@onid.orst.edu mailto:musta...@onid.orst.edu wrote: Dawei, I have no problems with proteins in the thousands of atoms. Can you post your command line and mdp files? Thank you, TJ Mustard On February 7, 2011 at 9:31 AM Da-Wei Li lida...@gmail.com mailto:lida...@gmail.com wrote: Well. It actually isn't dead but becomes very slow for large proteins. dawei On Mon, Feb 7, 2011 at 11:44 AM, Da-Wei Li lida...@gmail.com mailto:lida...@gmail.com wrote: hi, I did more test and found that it depended on size of the protein. Grompp will die when number of atoms of the protein is larger than about 200. Is it possible the source code limit the size of the protein that can be decoupled? thanks. dawei On Mon, Feb 7, 2011 at 10:34 AM, Da-Wei Li lida...@gmail.com mailto:lida...@gmail.com wrote: Dear users I tried free energy calculation but grompp couldn't go through. It stops after *** Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein' turning H bonds into constraints... Excluding 2 bonded neighbours molecule type 'SOL' turning H bonds into constraints...
Re: [gmx-users] Re: free energy calculation , grompp crash
Da-Wei, Do you need FEP information on the PR step? Are you going to do a MD(sd) with FEP on after the PR? And are you doing hydration of a protein? Thank you, TJ Mustard On February 7, 2011 at 10:23 AM Justin A. Lemkul jalem...@vt.edu wrote: Da-Wei Li wrote: hello Here they are the command line and mdp file. I use Gromacs 4.5.3. This is a test case only and the protein is 1UBQ. Grompp wills top for about 10 minutes then go through. The efficiency of this kind of process will depend on the amount of available memory on the system. Youre asking grompp to decouple a huge amount of degrees of freedom, which will require a lot of memory to do. It sounds like its working, in any case, so theres no real problem. Whether or not simultaneously decoupling the LJ and Coulombic interactions of a whole protein will generate a stable trajectory or sensible result is another matter. -Justin ***output of grompp* Back Off! I just backed up mdout.mdp to ./#mdout.mdp.6# Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type Protein turning H bonds into constraints... Excluding 2 bonded neighbours molecule type SOL turning H bonds into constraints... Coupling 1 copies of molecule type Protein Setting gen_seed to 8552 Velocities were taken from a Maxwell distribution at 300 K Command line and mdp file: ** grompp -f pr1.mdp -c after_em.gro -t em.trr -p topol.top -o pr1 ** define = -DPOSRES ; position restrain the protein ; Run parameters integrator = sd ; leap-frog integrator nsteps = 5000 ; 2 * 5 = 100 ps dt = 0.002 ; 2 fs ; Output control nstxout = 1000 ; save coordinates every 2 ps nstvout = 5000 ; save velocities every 100ps nstenergy = 1000 ; save energies every 2 ps nstlog = 1000 ; update log file every 2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = hbonds ; H bonds constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cels nstlist = 10 ; 20 fs rlist = 0.8 ; short-range neighborlist cutoff (in nm) rcoulomb = 0.8 ; short-range electrostatic cutoff (in nm) rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein Non-Protein ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no ; no pressure coupling in NVT ; 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 = 300 ; temperature for Maxwell distribution gen_seed = -1 ; generate a random seed ;free energy stuff free_energy = yes init_lambda = 0.0 delta_lambda = 0 sc_alpha=0.5 sc-power=1.0 sc-sigma= 0.3 couple-moltype = Protein couple-lambda0 = vdw-q couple-lambda1 = none *** thanks. dawei On Mon, Feb 7, 2011 at 1:05 PM, TJ Mustard musta...@onid.orst.edu mailto:musta...@onid.orst.edu wrote: Dawei, I have no problems with proteins in the thousands of atoms. Can you post your command line and mdp files? Thank you, TJ Mustard On February 7, 2011 at 9:31 AM Da-Wei Li lida...@gmail.com mailto:lida...@gmail.com wrote: Well. It actually isnt dead but becomes very slow for large
Re: [gmx-users] Re: free energy calculation , grompp crash
hello My goal is to study the solvation free energy of a fixed protein and compare it with implicit model. The pr1. mdp is just a test case. Grompp always need more than 10 minutes to finish for my 76 residues protein when free_energy = yes, no matter there is PR or not, whether I switch off only vdm or switch both ele and vdm. The memory usage is 1% on a system with 16GB memory so that memory limitation can be ruled out. best, dawei On Mon, Feb 7, 2011 at 1:31 PM, TJ Mustard musta...@onid.orst.edu wrote: Da-Wei, Do you need FEP information on the PR step? Are you going to do a MD(sd) with FEP on after the PR? And are you doing hydration of a protein? Thank you, TJ Mustard On February 7, 2011 at 10:23 AM Justin A. Lemkul jalem...@vt.edu wrote: Da-Wei Li wrote: hello Here they are the command line and mdp file. I use Gromacs 4.5.3. This is a test case only and the protein is 1UBQ. Grompp wills top for about 10 minutes then go through. The efficiency of this kind of process will depend on the amount of available memory on the system. You're asking grompp to decouple a huge amount of degrees of freedom, which will require a lot of memory to do. It sounds like it's working, in any case, so there's no real problem. Whether or not simultaneously decoupling the LJ and Coulombic interactions of a whole protein will generate a stable trajectory or sensible result is another matter. -Justin ***output of grompp* Back Off! I just backed up mdout.mdp to ./#mdout.mdp.6# Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein' turning H bonds into constraints... Excluding 2 bonded neighbours molecule type 'SOL' turning H bonds into constraints... Coupling 1 copies of molecule type 'Protein' Setting gen_seed to 8552 Velocities were taken from a Maxwell distribution at 300 K Command line and mdp file: ** grompp -f pr1.mdp -c after_em.gro -t em.trr -p topol.top -o pr1 ** define = -DPOSRES ; position restrain the protein ; Run parameters integrator = sd ; leap-frog integrator nsteps = 5000 ; 2 * 5 = 100 ps dt = 0.002 ; 2 fs ; Output control nstxout = 1000 ; save coordinates every 2 ps nstvout = 5000 ; save velocities every 100ps nstenergy = 1000 ; save energies every 2 ps nstlog = 1000 ; update log file every 2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = hbonds ; H bonds constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cels nstlist = 10 ; 20 fs rlist = 0.8 ; short-range neighborlist cutoff (in nm) rcoulomb = 0.8 ; short-range electrostatic cutoff (in nm) rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein Non-Protein ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no ; no pressure coupling in NVT ; 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 = 300 ; temperature for Maxwell distribution gen_seed = -1 ; generate a random seed ;free energy stuff free_energy = yes init_lambda = 0.0 delta_lambda = 0 sc_alpha =0.5 sc-power =1.0 sc-sigma = 0.3 couple-moltype = Protein couple-lambda0 = vdw-q couple-lambda1 = none *** thanks. dawei On Mon, Feb 7, 2011 at 1:05 PM, TJ Mustard musta...@onid.orst.edu mailto:musta...@onid.orst.edu wrote: Dawei, I have no problems with proteins in the thousands of atoms. Can you post your command line and mdp files? Thank you, TJ Mustard On February 7, 2011 at 9:31 AM Da-Wei Li lida...@gmail.com mailto:lida...@gmail.com wrote: Well. It actually isn't dead but becomes very slow for large proteins. dawei On Mon, Feb 7, 2011 at
Re: [gmx-users] Re: free energy calculation , grompp crash
hi, all I get another strange error. Once I turn on free energy calculation. mdrun will complain this: Initial maximum inter charge-group distances: two-body bonded interactions: 4.509 nm, LJC Pairs NB, atoms 1013 1231 multi-body bonded interactions: 0.428 nm, Proper Dih., atoms 5 13 Minimum cell size due to bonded interactions: 4.960 nm so that domain decomposition will not work. What is a LJC pair? I did not see any problem with my system. Everything is fine if free_energy = no. thanks. dawei On Mon, Feb 7, 2011 at 1:38 PM, Da-Wei Li lida...@gmail.com wrote: hello My goal is to study the solvation free energy of a fixed protein and compare it with implicit model. The pr1. mdp is just a test case. Grompp always need more than 10 minutes to finish for my 76 residues protein when free_energy = yes, no matter there is PR or not, whether I switch off only vdm or switch both ele and vdm. The memory usage is 1% on a system with 16GB memory so that memory limitation can be ruled out. best, dawei On Mon, Feb 7, 2011 at 1:31 PM, TJ Mustard musta...@onid.orst.edu wrote: Da-Wei, Do you need FEP information on the PR step? Are you going to do a MD(sd) with FEP on after the PR? And are you doing hydration of a protein? Thank you, TJ Mustard On February 7, 2011 at 10:23 AM Justin A. Lemkul jalem...@vt.edu wrote: Da-Wei Li wrote: hello Here they are the command line and mdp file. I use Gromacs 4.5.3. This is a test case only and the protein is 1UBQ. Grompp wills top for about 10 minutes then go through. The efficiency of this kind of process will depend on the amount of available memory on the system. You're asking grompp to decouple a huge amount of degrees of freedom, which will require a lot of memory to do. It sounds like it's working, in any case, so there's no real problem. Whether or not simultaneously decoupling the LJ and Coulombic interactions of a whole protein will generate a stable trajectory or sensible result is another matter. -Justin ***output of grompp* Back Off! I just backed up mdout.mdp to ./#mdout.mdp.6# Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein' turning H bonds into constraints... Excluding 2 bonded neighbours molecule type 'SOL' turning H bonds into constraints... Coupling 1 copies of molecule type 'Protein' Setting gen_seed to 8552 Velocities were taken from a Maxwell distribution at 300 K Command line and mdp file: ** grompp -f pr1.mdp -c after_em.gro -t em.trr -p topol.top -o pr1 ** define = -DPOSRES ; position restrain the protein ; Run parameters integrator = sd ; leap-frog integrator nsteps = 5000 ; 2 * 5 = 100 ps dt = 0.002 ; 2 fs ; Output control nstxout = 1000 ; save coordinates every 2 ps nstvout = 5000 ; save velocities every 100ps nstenergy = 1000 ; save energies every 2 ps nstlog = 1000 ; update log file every 2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = hbonds ; H bonds constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cels nstlist = 10 ; 20 fs rlist = 0.8 ; short-range neighborlist cutoff (in nm) rcoulomb = 0.8 ; short-range electrostatic cutoff (in nm) rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein Non-Protein ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no ; no pressure coupling in NVT ; 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 = 300 ; temperature for Maxwell distribution gen_seed = -1 ; generate a random seed ;free energy stuff free_energy = yes init_lambda = 0.0 delta_lambda = 0 sc_alpha =0.5 sc-power =1.0 sc-sigma = 0.3 couple-moltype = Protein couple-lambda0 = vdw-q couple-lambda1 = none
Re: [gmx-users] Re: free energy calculation , grompp crash
Da-Wei, How do you generate your box, solvent, ions? I use this sequence of commands: pdb2gmx -f protein.pdb -o protein.gro -p protein.top editconf -bt cubic -f protein.gro -o protein.gro -c -d 1.5 genbox -cp protein.gro -cs spc216.gro -o protein_b4ion.gro -p protein.top grompp -f em.mdp -c protein_b4ion.gro -p protein.top -o protein_b4ion.tpr genion -s protein_b4ion.tpr -o protein_b4em.gro -neutral -conc 0.001 -pname NA -nname CL -g protein_ion.log -p protein.top Then I will grompp and mdrun. I seldom change these values. Thank you, TJ Musard On February 7, 2011 at 1:08 PM Da-Wei Li lida...@gmail.com wrote: hi, all I get another strange error. Once I turn on free energy calculation. mdrun will complain this: Initial maximum inter charge-group distances: two-body bonded interactions: 4.509 nm, LJC Pairs NB, atoms 1013 1231 multi-body bonded interactions: 0.428 nm, Proper Dih., atoms 5 13 Minimum cell size due to bonded interactions: 4.960 nm so that domain decomposition will not work. What is a LJC pair? I did not see any problem with my system. Everything is fine iffree_energy = no. thanks. dawei On Mon, Feb 7, 2011 at 1:38 PM, Da-Wei Li lida...@gmail.com wrote: hello My goal is to study the solvation free energy of a fixed protein and compare it with implicit model. The pr1. mdp is just a test case. Grompp always need more than 10 minutes to finish for my 76 residues protein when free_energy = yes, no matter there is PR or not, whether I switch off only vdm or switch both ele and vdm. The memory usage is 1% on a system with 16GB memory so that memory limitation can be ruled out. best, dawei On Mon, Feb 7, 2011 at 1:31 PM, TJ Mustard musta...@onid.orst.edu wrote: Da-Wei, Do you need FEP information on the PR step? Are you going to do a MD(sd) with FEP on after the PR? And are you doing hydration of a protein? Thank you, TJ Mustard On February 7, 2011 at 10:23 AM Justin A. Lemkul jalem...@vt.edu wrote: Da-Wei Li wrote: hello Here they are the command line and mdp file. I use Gromacs 4.5.3. This is a test case only and the protein is 1UBQ. Grompp wills top for about 10 minutes then go through. The efficiency of this kind of process will depend on the amount of available memory on the system. Youre asking grompp to decouple a huge amount of degrees of freedom, which will require a lot of memory to do. It sounds like its working, in any case, so theres no real problem. Whether or not simultaneously decoupling the LJ and Coulombic interactions of a whole protein will generate a stable trajectory or sensible result is another matter. -Justin ***output of grompp* Back Off! I just backed up mdout.mdp to ./#mdout.mdp.6# Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5
Re: [gmx-users] Re: free energy calculation , grompp crash
Da-Wei Li wrote: hi, all I get another strange error. Once I turn on free energy calculation. mdrun will complain this: Initial maximum inter charge-group distances: two-body bonded interactions: 4.509 nm, LJC Pairs NB, atoms 1013 1231 multi-body bonded interactions: 0.428 nm, Proper Dih., atoms 5 13 Minimum cell size due to bonded interactions: 4.960 nm It looks like your topology has some serious problems. Bonded interactions should not occur at 4-nm distance, unless you're using some special distance restraints, in which case you probably can't make DD work very efficiently; mdrun -pd might be an option in that case. If you're running a basic protein simulation without special restraints, etc then something is wrong in the topology. so that domain decomposition will not work. What is a LJC pair? I did Some sort of 1-4 interaction, it looks like, but I didn't spend much time in the code to be absolutely sure. -Justin not see any problem with my system. Everything is fine if free_energy = no. thanks. dawei On Mon, Feb 7, 2011 at 1:38 PM, Da-Wei Li lida...@gmail.com mailto:lida...@gmail.com wrote: hello My goal is to study the solvation free energy of a fixed protein and compare it with implicit model. The pr1. mdp is just a test case. Grompp always need more than 10 minutes to finish for my 76 residues protein when free_energy = yes, no matter there is PR or not, whether I switch off only vdm or switch both ele and vdm. The memory usage is 1% on a system with 16GB memory so that memory limitation can be ruled out. best, dawei On Mon, Feb 7, 2011 at 1:31 PM, TJ Mustard musta...@onid.orst.edu mailto:musta...@onid.orst.edu wrote: Da-Wei, Do you need FEP information on the PR step? Are you going to do a MD(sd) with FEP on after the PR? And are you doing hydration of a protein? Thank you, TJ Mustard On February 7, 2011 at 10:23 AM Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: Da-Wei Li wrote: hello Here they are the command line and mdp file. I use Gromacs 4.5.3. This is a test case only and the protein is 1UBQ. Grompp wills top for about 10 minutes then go through. The efficiency of this kind of process will depend on the amount of available memory on the system. You're asking grompp to decouple a huge amount of degrees of freedom, which will require a lot of memory to do. It sounds like it's working, in any case, so there's no real problem. Whether or not simultaneously decoupling the LJ and Coulombic interactions of a whole protein will generate a stable trajectory or sensible result is another matter. -Justin ***output of grompp* Back Off! I just backed up mdout.mdp to ./#mdout.mdp.6# Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein' turning H bonds into constraints... Excluding 2 bonded neighbours molecule type 'SOL' turning H bonds into constraints... Coupling 1 copies of molecule type 'Protein' Setting gen_seed to 8552 Velocities were taken from a Maxwell distribution at 300 K Command line and mdp file: ** grompp -f pr1.mdp -c after_em.gro -t em.trr -p topol.top -o pr1 ** define = -DPOSRES ; position restrain the protein ; Run parameters integrator = sd ; leap-frog integrator nsteps = 5000 ; 2 * 5 = 100 ps dt = 0.002 ; 2 fs ; Output control nstxout = 1000 ; save coordinates every 2 ps nstvout = 5000 ; save velocities every 100ps nstenergy = 1000 ; save energies every 2 ps nstlog = 1000 ; update log file every 2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = hbonds ; H bonds constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cels nstlist = 10 ; 20 fs rlist = 0.8 ; short-range neighborlist cutoff (in nm) rcoulomb = 0.8 ; short-range
Re: [gmx-users] Re: free energy calculation , grompp crash
hello 1. I used almost identical comands with you. 2. If I set couple-intramol = yes. My grompp will run very fast. From the manu, couple-intramol = no means all intramol (protein-protein interaction in my case) interaction are turned to predefind list while couple-intramol = yes means intra-protein interactions are scaled the same with protein-water interaction. Maybe this is the reason my grompp becomes very slow?? 3. If I set couple-intramol = yes, my mdrun will go very smooth. It is also fine if I turn off free energy calculation. If I set couple-intramol=no, it will complain two-body bonded interactions: 4.509 nm, LJC Pairs NB, atoms 1013 1231. This is very strange. Atom 1013 and 1231 are not bonded at all, they are actually the most *far away* pair in my system. I can't understand why free energy calculation will cause this message. best, dawei On Mon, Feb 7, 2011 at 6:16 PM, TJ Mustard musta...@onid.orst.edu wrote: Da-Wei, How do you generate your box, solvent, ions? I use this sequence of commands: pdb2gmx -f protein.pdb -o protein.gro -p protein.top editconf -bt cubic -f protein.gro -o protein.gro -c -d 1.5 genbox -cp protein.gro -cs spc216.gro -o protein_b4ion.gro -p protein.top grompp -f em.mdp -c protein_b4ion.gro -p protein.top -o protein_b4ion.tpr genion -s protein_b4ion.tpr -o protein_b4em.gro -neutral -conc 0.001 -pname NA -nname CL -g protein_ion.log -p protein.top Then I will grompp and mdrun. I seldom change these values. Thank you, TJ Musard On February 7, 2011 at 1:08 PM Da-Wei Li lida...@gmail.com wrote: hi, all I get another strange error. Once I turn on free energy calculation. mdrun will complain this: Initial maximum inter charge-group distances: two-body bonded interactions: 4.509 nm, LJC Pairs NB, atoms 1013 1231 multi-body bonded interactions: 0.428 nm, Proper Dih., atoms 5 13 Minimum cell size due to bonded interactions: 4.960 nm so that domain decomposition will not work. What is a LJC pair? I did not see any problem with my system. Everything is fine if free_energy = no. thanks. dawei On Mon, Feb 7, 2011 at 1:38 PM, Da-Wei Li lida...@gmail.com wrote: hello My goal is to study the solvation free energy of a fixed protein and compare it with implicit model. The pr1. mdp is just a test case. Grompp always need more than 10 minutes to finish for my 76 residues protein when free_energy = yes, no matter there is PR or not, whether I switch off only vdm or switch both ele and vdm. The memory usage is 1% on a system with 16GB memory so that memory limitation can be ruled out. best, dawei On Mon, Feb 7, 2011 at 1:31 PM, TJ Mustard musta...@onid.orst.edu wrote: Da-Wei, Do you need FEP information on the PR step? Are you going to do a MD(sd) with FEP on after the PR? And are you doing hydration of a protein? Thank you, TJ Mustard On February 7, 2011 at 10:23 AM Justin A. Lemkul jalem...@vt.edu wrote: Da-Wei Li wrote: hello Here they are the command line and mdp file. I use Gromacs 4.5.3. This is a test case only and the protein is 1UBQ. Grompp wills top for about 10 minutes then go through. The efficiency of this kind of process will depend on the amount of available memory on the system. You're asking grompp to decouple a huge amount of degrees of freedom, which will require a lot of memory to do. It sounds like it's working, in any case, so there's no real problem. Whether or not simultaneously decoupling the LJ and Coulombic interactions of a whole protein will generate a stable trajectory or sensible result is another matter. -Justin ***output of grompp* Back Off! I just backed up mdout.mdp to ./#mdout.mdp.6# Generated 2278 of the 2278 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 2278 of the 2278 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein' turning H bonds into constraints... Excluding 2 bonded neighbours molecule type 'SOL' turning H bonds into constraints... Coupling 1 copies of molecule type 'Protein' Setting gen_seed to 8552 Velocities were taken from a Maxwell distribution at 300 K Command line and mdp file: ** grompp -f pr1.mdp -c after_em.gro -t em.trr -p topol.top -o pr1 ** define = -DPOSRES ; position restrain the protein ; Run parameters integrator = sd ; leap-frog integrator nsteps = 5000 ; 2 * 5 = 100 ps dt = 0.002 ; 2 fs ; Output control nstxout = 1000 ; save coordinates every 2 ps nstvout = 5000 ; save velocities every 100ps nstenergy = 1000 ; save energies every 2 ps nstlog = 1000 ; update log file every 2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic
Re: [gmx-users] Re: free energy calculation , grompp crash
Da-Wei Li wrote: hello 1. I used almost identical comands with you. 2. If I set couple-intramol = yes. My grompp will run very fast. From the manu, couple-intramol = no means all intramol (protein-protein interaction in my case) interaction are turned to predefind list while couple-intramol = yes means intra-protein interactions are scaled the same with protein-water interaction. Maybe this is the reason my grompp becomes very slow?? That seems likely. As I suggested before, grompp is chewing up a lot of memory converting all intramolecular nonbonded interactions into pairs and exclusions. This is a very intensive process, and scales with system size. 3. If I set couple-intramol = yes, my mdrun will go very smooth. It is also fine if I turn off free energy calculation. If I set couple-intramol=no, it will complain two-body bonded interactions: 4.509 nm, LJC Pairs NB, atoms 1013 1231. This is very strange. Atom 1013 and 1231 are not bonded at all, they are actually the most *far away* pair in my system. I can't understand why free energy calculation will cause this message. This is a consequence of the description above. Check the output of grompp -pp or the .tpr file and you will likely see that interactions within your protein are now being assigned to the [pairs] and [exclusions] directives. As a result, you have bonded interactions that are very far apart. You might be able to use mdrun -pd to circumvent this performance issue (although it will be slower than DD in most cases, here it will likely be faster than any construct DD might come up with, if it runs at all). -Justin best, dawei On Mon, Feb 7, 2011 at 6:16 PM, TJ Mustard musta...@onid.orst.edu mailto:musta...@onid.orst.edu wrote: Da-Wei, How do you generate your box, solvent, ions? I use this sequence of commands: pdb2gmx -f protein.pdb -o protein.gro -p protein.top editconf -bt cubic -f protein.gro -o protein.gro -c -d 1.5 genbox -cp protein.gro -cs spc216.gro -o protein_b4ion.gro -p protein.top grompp -f em.mdp -c protein_b4ion.gro -p protein.top -o protein_b4ion.tpr genion -s protein_b4ion.tpr -o protein_b4em.gro -neutral -conc 0.001 -pname NA -nname CL -g protein_ion.log -p protein.top Then I will grompp and mdrun. I seldom change these values. Thank you, TJ Musard On February 7, 2011 at 1:08 PM Da-Wei Li lida...@gmail.com mailto:lida...@gmail.com wrote: hi, all I get another strange error. Once I turn on free energy calculation. mdrun will complain this: Initial maximum inter charge-group distances: two-body bonded interactions: 4.509 nm, LJC Pairs NB, atoms 1013 1231 multi-body bonded interactions: 0.428 nm, Proper Dih., atoms 5 13 Minimum cell size due to bonded interactions: 4.960 nm so that domain decomposition will not work. What is a LJC pair? I did not see any problem with my system. Everything is fine if free_energy = no. thanks. dawei On Mon, Feb 7, 2011 at 1:38 PM, Da-Wei Li lida...@gmail.com mailto:lida...@gmail.com wrote: hello My goal is to study the solvation free energy of a fixed protein and compare it with implicit model. The pr1. mdp is just a test case. Grompp always need more than 10 minutes to finish for my 76 residues protein when free_energy = yes, no matter there is PR or not, whether I switch off only vdm or switch both ele and vdm. The memory usage is 1% on a system with 16GB memory so that memory limitation can be ruled out. best, dawei On Mon, Feb 7, 2011 at 1:31 PM, TJ Mustard musta...@onid.orst.edu mailto:musta...@onid.orst.edu wrote: Da-Wei, Do you need FEP information on the PR step? Are you going to do a MD(sd) with FEP on after the PR? And are you doing hydration of a protein? Thank you, TJ Mustard On February 7, 2011 at 10:23 AM Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: Da-Wei Li wrote: hello Here they are the command line and mdp file. I use Gromacs 4.5.3. This is a test case only and the protein is 1UBQ. Grompp wills top for about 10 minutes then go through. The efficiency of this kind of process will depend on the amount of available memory on the system. You're asking grompp to decouple a huge amount of degrees of freedom, which will require a lot of memory to do. It sounds like it's working, in any case, so there's no real problem.
[gmx-users] Re: Free energy calculations - desolvation energy of Na+
Hi Justin, thanks for this link: http://lists.gromacs.org/pipermail/gmx-users/2003-June/006054.html It clearly says, why sampling by MD is not proper. Best, Eva Hello all, I am trying calculate desolvation free energy of Na+ in water using option couple-moltype, not by creating B-topology in the .top file. At first, I switched off coulombic interaction using: couple-moltype = Na+ couple-lambda0 = vdw-q couple-lambda1 = vdw Then I tried to switch off LJ interaction using couple-lambda0 = vdw and couple-lambda1 = none, but obtained this warning: WARNING 1 [file equil_NVT000.mdp, line 95]: For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used This is strange, because my .mdp file has only 94 lines. Regardless of the line number, the error message is pretty explicit - you probably have integrator = md in your .mdp file when you should have integrator = sd. -Justin Thanks in advance for any help or suggestion. Best, Eva For convenience, I am pasting weak coupling algorithm section and free energy section from my .mdp file. The simulation was performed at constant volume using strong temperature coupling, since it was an equilibration run. ; OPTIONS FOR WEAK COUPLING ALGORITHMS Tcoupl = v-rescale tc-grps = system tau_t= 0.1 ref_t= 300 nsttcouple = 1 ; Free energy control stuff free_energy = yes init_lambda = 0.0 delta_lambda = 0 foreign_lambda = 0.1 sc-alpha = 0.7 sc-power = 1 sc-sigma = 0.3 nstdhdl = 10 separate-dhdl-file = yes dhdl-derivatives = yes dh_hist_size = 0 dh_hist_spacing = 0.1 couple-moltype = Na+ couple-lambda0 = vdw couple-lambda1 = none couple-intramol = no -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- -- gmx-users mailing list gmx-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! End of gmx-users Digest, Vol 78, Issue 186 ** -- 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] Re: Free Energy Calculation
The question is: Do you trust you simulations when you/others do not find significant problems or do you only trust your simulations when you are sure the integrators etc are implemented exactly correctly or at least the magnitude of the (small) errors are known. Personally I (and the majority of people) fall in the first category. David is clearly in the second category. I think you can recognise free energy calculation specialists like David and myself by a heightened level of paranoia, along with other signs like staring blankly at the wall wondering which part of your protocol could be causing that 5 kcal/mol deviation... so while I'm a big supporter of the first category I see where David's coming from, having spent a lot of time in the past running controls for almost certainly inconsequential variables like this. If I can offer my two eurocents-worth, I would sleep easier at night with a barostat that's both well-behaved and produces a rigorous NPT ensemble... I'm in the process of migrating to GROMACS from NAMD and have never heard any complaints about their Langevin piston thermostat (http://dx.doi.org/10.1063/1.470648), would there be any downsides to having this method in GROMACS? best wishes, Floris From: Berk Hess g...@hotmail.com To: Discussion list for GROMACS users gmx-users@gromacs.org Sent: Friday, 7 August, 2009 16:16:19 Subject: RE: [gmx-users] Re: Free Energy Calculation All the discussion is mainly a matter about how paranoid you are. A Berendsen barostat does not generate a proper ensemble. But it does get the average pressure right. The questions is how large the effects of the ensemble error is. I would say that in practice this would be negligible compared to fluctuations, force field errors, insufficient sampling etc. But strictly speaking it is there. Parrinello-Rahman in principle gives you a proper NPT ensemble (combined with a proper thermostat). But to have a real symplectic or even only reversible integrator is difficult and might be computationally expensive, especially when constraints are present, as will nearly always be the case with most water models. Again errors of a non-reversible integrator are, in most cases, extremely small and smaller than the other errors you make. In addition PR-coupling might produce nasty volume oscillations. The question is: Do you trust you simulations when you/others do not find significant problems or do you only trust your simulations when you are sure the integrators etc are implemented exactly correctly or at least the magnitude of the (small) errors are known. Personally I (and the majority of people) fall in the first category. David is clearly in the second category. In the past decade it has been shown that several algorithms that had been used a lot caused serious artifacts. But I would say that we are currently in a state where most things have been checked and although we know there are some errors (BTW descretizing analytic differential equations also produces errors), these errors have no effect on MOST of the things we are interested in. Of course you should not derive compressibilities using a Berendsen barostat, as a big name in the field has done even not to long ago. Berk Date: Fri, 7 Aug 2009 07:32:33 -0700 From: floris_buel...@yahoo.com Subject: Re: [gmx-users] Re: Free Energy Calculation To: gmx-users@gromacs.org Hi David and Berk, I guess the best proper solution would be to use a Parrinello-Rahman barostat. But in Gromacs this is currently not implemented in a reversible way. Yes, I would love to see a barostat that samples the correct distribution. Can you clarify? The wording in the manual implies that Parrinello-Rahman coupling samples the correct distribution... When you say it's not reversible do I understand right that means you can't run your simulation backwards and get an identical result, if yes, is this ever relevant in practice for e.g. solvation free energy calculation? David, why don't you use Parrinello-Rahman coupling? thanks, Floris In practice I think the issues with the Berendsen thermostat are not relevant. It does not generate the correct pressure and volume fluctuations, but the average pressure and volume are correct. The effects of incorrect fluctuations is negligible and is surely much smaller than simulating at the wrong volume, which can have a large effect on solvation free energies. It seems fairly like you are right (assuming your first line is meant to refer to the Berendsen barostat), but I think no one knows for sure at this point whether it matters and when it might matter. I agree that simulating at an incorrect volume can be a big problem, however. This is unlikely, but can happen, which is why we now do the affine transformation I just mentioned when looking at solvation free energies. Anyway -- I'd just like to see a proper barostat implemented
RE: [gmx-users] Re: Free Energy Calculation
Good question. I guess the best proper solution would be to use a Parrinello-Rahman barostat. But in Gromacs this is currently not implemented in a reversible way. In practice I think the issues with the Berendsen thermostat are not relevant. It does not generate the correct pressure and volume fluctuations, but the average pressure and volume are correct. The effects of incorrect fluctuations is negligible and is surely much smaller than simulating at the wrong volume, which can have a large effect on solvation free energies. I also think that issues with Berendsen pressure coupling are much less severe than the issues with Berendsen temperature coupling. Berendsen tcoupl can have significant effects on your ensemble. But in Gromacs 4.0 this can easily be avoided by using the v-rescale thermostat. Berk Subject: RE: [gmx-users] Re: Free Energy Calculation Date: Thu, 6 Aug 2009 16:01:04 -0400 From: pn...@utnet.utoledo.edu To: gmx-users@gromacs.org RE: [gmx-users] Re: Free Energy Calculation A question to the answer below: If the Berendsen barostat does not sample correctly the distribution of pressures then what guarantees that the the volume/density is correct at the end of the equilibration phase when the user is advised to switch to NVT? Peter Nagy The University of Toledo -Original Message- From: gmx-users-boun...@gromacs.org on behalf of David Mobley Sent: Thu 8/6/2009 2:38 PM To: Mauricio Carrillo Tripp Cc: Discussion list for GROMACS users Subject: [gmx-users] Re: Free Energy Calculation Hi, Sorry for the delay answering. These questions are better put on the GROMACS users list. 1) Yes -- the Berendsen barostat does not sample the correct distribution of pressures. 2) Regenerating velocities is fine for a couple of reasons: (a) we are after thermodynamics, not dynamics, and (b) regenerating velocities does not de-equilibrate the system (see the Andersen thermostat, for example). On Jun 8, 2009, at 3:02 PM, Mauricio Carrillo Tripp wrote: Hi, I've been following directions at http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial to calculate hydration free energies of a few small molecules, and I have a couple of questions: 1) Is there a reason for using NVT instead of NPT for the production step? 2) I noticed that between the equilibration and the production runs you generate velocities at each step, instead of using the ones from the previous step. Is there a reason for doing this? wouldn't this 'de-equilibrate' an already equilibrated system? Thank you. -- 0 | Mauricio Carrillo Tripp, PhD / | Department of Molecular Biology, TPC6 0| The Scripps Research Institute \ | 10550 North Torrey Pines Road 0 | La Jolla, California 92037 / | tri...@scripps.edu 0| http://www.scripps.edu/~trippm ** Aut tace aut loquere meliora silentio ** _ What can you do with the new Windows Live? Find out http://www.microsoft.com/windows/windowslive/default.aspx___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Re: Free Energy Calculation
Peter, On Thu, Aug 6, 2009 at 3:01 PM, Nagy, Peter I. pn...@utnet.utoledo.eduwrote: A question to the answer below: If the Berendsen barostat does not sample correctly the distribution of pressures then what guarantees that the the volume/density is correct at the end of the equilibration phase when the user is advised to switch to NVT? Very good question. In fact, nothing. More recently, we do an affine transformation to rescale the box size to the correct (average) box size at the end of NPT equilibration and then run a short additional NVT eqiulibration before NVT production. However, this tends not to work well for systems containing large molecules such as proteins. Ideally I would instead have a barostat that samples the proper distribution of pressures to use! David Peter Nagy The University of Toledo -Original Message- From: gmx-users-boun...@gromacs.org on behalf of David Mobley Sent: Thu 8/6/2009 2:38 PM To: Mauricio Carrillo Tripp Cc: Discussion list for GROMACS users Subject: [gmx-users] Re: Free Energy Calculation Hi, Sorry for the delay answering. These questions are better put on the GROMACS users list. 1) Yes -- the Berendsen barostat does not sample the correct distribution of pressures. 2) Regenerating velocities is fine for a couple of reasons: (a) we are after thermodynamics, not dynamics, and (b) regenerating velocities does not de-equilibrate the system (see the Andersen thermostat, for example). On Jun 8, 2009, at 3:02 PM, Mauricio Carrillo Tripp wrote: Hi, I've been following directions at http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial to calculate hydration free energies of a few small molecules, and I have a couple of questions: 1) Is there a reason for using NVT instead of NPT for the production step? 2) I noticed that between the equilibration and the production runs you generate velocities at each step, instead of using the ones from the previous step. Is there a reason for doing this? wouldn't this 'de-equilibrate' an already equilibrated system? Thank you. -- 0 | Mauricio Carrillo Tripp, PhD / | Department of Molecular Biology, TPC6 0| The Scripps Research Institute \ | 10550 North Torrey Pines Road 0 | La Jolla, California 92037 / | tri...@scripps.edu 0| http://www.scripps.edu/~trippmhttp://www.scripps.edu/%7Etrippm ** Aut tace aut loquere meliora silentio ** ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Re: Free Energy Calculation
Berk, I guess the best proper solution would be to use a Parrinello-Rahman barostat. But in Gromacs this is currently not implemented in a reversible way. Yes, I would love to see a barostat that samples the correct distribution. In practice I think the issues with the Berendsen thermostat are not relevant. It does not generate the correct pressure and volume fluctuations, but the average pressure and volume are correct. The effects of incorrect fluctuations is negligible and is surely much smaller than simulating at the wrong volume, which can have a large effect on solvation free energies. It seems fairly like you are right (assuming your first line is meant to refer to the Berendsen barostat), but I think no one knows for sure at this point whether it matters and when it might matter. I agree that simulating at an incorrect volume can be a big problem, however. This is unlikely, but can happen, which is why we now do the affine transformation I just mentioned when looking at solvation free energies. Anyway -- I'd just like to see a proper barostat implemented at some point, then the paranoid (like me) can run with a proper barostat instead of worrying about whether or not this matters, AND someone could test whether observables they are interested in come out different depending on the choice of barostat. David I also think that issues with Berendsen pressure coupling are much less severe than the issues with Berendsen temperature coupling. Berendsen tcoupl can have significant effects on your ensemble. But in Gromacs 4.0 this can easily be avoided by using the v-rescale thermostat. Berk -- Subject: RE: [gmx-users] Re: Free Energy Calculation Date: Thu, 6 Aug 2009 16:01:04 -0400 From: pn...@utnet.utoledo.edu To: gmx-users@gromacs.org A question to the answer below: If the Berendsen barostat does not sample correctly the distribution of pressures then what guarantees that the the volume/density is correct at the end of the equilibration phase when the user is advised to switch to NVT? Peter Nagy The University of Toledo -Original Message- From: gmx-users-boun...@gromacs.org on behalf of David Mobley Sent: Thu 8/6/2009 2:38 PM To: Mauricio Carrillo Tripp Cc: Discussion list for GROMACS users Subject: [gmx-users] Re: Free Energy Calculation Hi, Sorry for the delay answering. These questions are better put on the GROMACS users list. 1) Yes -- the Berendsen barostat does not sample the correct distribution of pressures. 2) Regenerating velocities is fine for a couple of reasons: (a) we are after thermodynamics, not dynamics, and (b) regenerating velocities does not de-equilibrate the system (see the Andersen thermostat, for example). On Jun 8, 2009, at 3:02 PM, Mauricio Carrillo Tripp wrote: Hi, I've been following directions at http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial to calculate hydration free energies of a few small molecules, and I have a couple of questions: 1) Is there a reason for using NVT instead of NPT for the production step? 2) I noticed that between the equilibration and the production runs you generate velocities at each step, instead of using the ones from the previous step. Is there a reason for doing this? wouldn't this 'de-equilibrate' an already equilibrated system? Thank you. -- 0 | Mauricio Carrillo Tripp, PhD / | Department of Molecular Biology, TPC6 0| The Scripps Research Institute \ | 10550 North Torrey Pines Road 0 | La Jolla, California 92037 / | tri...@scripps.edu 0| http://www.scripps.edu/~trippmhttp://www.scripps.edu/%7Etrippm ** Aut tace aut loquere meliora silentio ** -- What can you do with the new Windows Live? Find outhttp://www.microsoft.com/windows/windowslive/default.aspx ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Re: Free Energy Calculation
Hi David and Berk, I guess the best proper solution would be to use a Parrinello-Rahman barostat. But in Gromacs this is currently not implemented in a reversible way. Yes, I would love to see a barostat that samples the correct distribution. Can you clarify? The wording in the manual implies that Parrinello-Rahman coupling samples the correct distribution... When you say it's not reversible do I understand right that means you can't run your simulation backwards and get an identical result, if yes, is this ever relevant in practice for e.g. solvation free energy calculation? David, why don't you use Parrinello-Rahman coupling? thanks, Floris In practice I think the issues with the Berendsen thermostat are not relevant. It does not generate the correct pressure and volume fluctuations, but the average pressure and volume are correct. The effects of incorrect fluctuations is negligible and is surely much smaller than simulating at the wrong volume, which can have a large effect on solvation free energies. It seems fairly like you are right (assuming your first line is meant to refer to the Berendsen barostat), but I think no one knows for sure at this point whether it matters and when it might matter. I agree that simulating at an incorrect volume can be a big problem, however. This is unlikely, but can happen, which is why we now do the affine transformation I just mentioned when looking at solvation free energies. Anyway -- I'd just like to see a proper barostat implemented at some point, then the paranoid (like me) can run with a proper barostat instead of worrying about whether or not this matters, AND someone could test whether observables they are interested in come out different depending on the choice of barostat. David I also think that issues with Berendsen pressure coupling are much less severe than the issues with Berendsen temperature coupling. Berendsen tcoupl can have significant effects on your ensemble. But in Gromacs 4.0 this can easily be avoided by using the v-rescale thermostat. Berk Subject: RE: [gmx-users] Re: Free Energy Calculation Date: Thu, 6 Aug 2009 16:01:04 -0400 From: pn...@utnet.utoledo.edu To: gmx-users@gromacs.org A question to the answer below: If the Berendsen barostat does not sample correctly the distribution of pressures then what guarantees that the the volume/density is correct at the end of the equilibration phase when the user is advised to switch to NVT? Peter Nagy The University of Toledo -Original Message- From: gmx-users-boun...@gromacs.org on behalf of David Mobley Sent: Thu 8/6/2009 2:38 PM To: Mauricio Carrillo Tripp Cc: Discussion list for GROMACS users Subject: [gmx-users] Re: Free Energy Calculation Hi, Sorry for the delay answering. These questions are better put on the GROMACS users list. 1) Yes -- the Berendsen barostat does not sample the correct distribution of pressures. 2) Regenerating velocities is fine for a couple of reasons: (a) we are after thermodynamics, not dynamics, and (b) regenerating velocities does not de-equilibrate the system (see the Andersen thermostat, for example). On Jun 8, 2009, at 3:02 PM, Mauricio Carrillo Tripp wrote: Hi, I've been following directions at http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial to calculate hydration free energies of a few small molecules, and I have a couple of questions: 1) Is there a reason for using NVT instead of NPT for the production step? 2) I noticed that between the equilibration and the production runs you generate velocities at each step, instead of using the ones from the previous step. Is there a reason for doing this? wouldn't this 'de-equilibrate' an already equilibrated system? Thank you. -- 0 | Mauricio Carrillo Tripp, PhD / | Department of Molecular Biology, TPC6 0| The Scripps Research Institute \ | 10550 North Torrey Pines Road 0 | La Jolla, California 92037 / | tri...@scripps.edu 0| http://www.scripps.edu/~trippm ** Aut tace aut loquere meliora silentio ** What can you do with the new Windows Live? Find out ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe
RE: [gmx-users] Re: Free Energy Calculation
All the discussion is mainly a matter about how paranoid you are. A Berendsen barostat does not generate a proper ensemble. But it does get the average pressure right. The questions is how large the effects of the ensemble error is. I would say that in practice this would be negligible compared to fluctuations, force field errors, insufficient sampling etc. But strictly speaking it is there. Parrinello-Rahman in principle gives you a proper NPT ensemble (combined with a proper thermostat). But to have a real symplectic or even only reversible integrator is difficult and might be computationally expensive, especially when constraints are present, as will nearly always be the case with most water models. Again errors of a non-reversible integrator are, in most cases, extremely small and smaller than the other errors you make. In addition PR-coupling might produce nasty volume oscillations. The question is: Do you trust you simulations when you/others do not find significant problems or do you only trust your simulations when you are sure the integrators etc are implemented exactly correctly or at least the magnitude of the (small) errors are known. Personally I (and the majority of people) fall in the first category. David is clearly in the second category. In the past decade it has been shown that several algorithms that had been used a lot caused serious artifacts. But I would say that we are currently in a state where most things have been checked and although we know there are some errors (BTW descretizing analytic differential equations also produces errors), these errors have no effect on MOST of the things we are interested in. Of course you should not derive compressibilities using a Berendsen barostat, as a big name in the field has done even not to long ago. Berk Date: Fri, 7 Aug 2009 07:32:33 -0700 From: floris_buel...@yahoo.com Subject: Re: [gmx-users] Re: Free Energy Calculation To: gmx-users@gromacs.org Hi David and Berk, I guess the best proper solution would be to use a Parrinello-Rahman barostat. But in Gromacs this is currently not implemented in a reversible way. Yes, I would love to see a barostat that samples the correct distribution. Can you clarify? The wording in the manual implies that Parrinello-Rahman coupling samples the correct distribution... When you say it's not reversible do I understand right that means you can't run your simulation backwards and get an identical result, if yes, is this ever relevant in practice for e.g. solvation free energy calculation? David, why don't you use Parrinello-Rahman coupling? thanks, Floris In practice I think the issues with the Berendsen thermostat are not relevant. It does not generate the correct pressure and volume fluctuations, but the average pressure and volume are correct. The effects of incorrect fluctuations is negligible and is surely much smaller than simulating at the wrong volume, which can have a large effect on solvation free energies. It seems fairly like you are right (assuming your first line is meant to refer to the Berendsen barostat), but I think no one knows for sure at this point whether it matters and when it might matter. I agree that simulating at an incorrect volume can be a big problem, however. This is unlikely, but can happen, which is why we now do the affine transformation I just mentioned when looking at solvation free energies. Anyway -- I'd just like to see a proper barostat implemented at some point, then the paranoid (like me) can run with a proper barostat instead of worrying about whether or not this matters, AND someone could test whether observables they are interested in come out different depending on the choice of barostat. David I also think that issues with Berendsen pressure coupling are much less severe than the issues with Berendsen temperature coupling. Berendsen tcoupl can have significant effects on your ensemble. But in Gromacs 4.0 this can easily be avoided by using the v-rescale thermostat. Berk Subject: RE: [gmx-users] Re: Free Energy Calculation Date: Thu, 6 Aug 2009 16:01:04 -0400 From: pn...@utnet.utoledo.edu To: gmx-users@gromacs.org A question to the answer below: If the Berendsen barostat does not sample correctly the distribution of pressures then what guarantees that the the volume/density is correct at the end of the equilibration phase when the user is advised to switch to NVT? Peter Nagy The University of Toledo -Original Message- From: gmx-users-boun...@gromacs.org on behalf of David Mobley Sent: Thu 8/6/2009 2:38 PM To: Mauricio Carrillo Tripp Cc: Discussion list for GROMACS users Subject: [gmx-users] Re: Free Energy Calculation Hi, Sorry for the delay answering. These questions are better put on the GROMACS users list. 1) Yes -- the Berendsen barostat does not sample the correct distribution of pressures. 2) Regenerating velocities is fine for a couple
[gmx-users] Re: Free Energy Calculation
Hi, Sorry for the delay answering. These questions are better put on the GROMACS users list. 1) Yes -- the Berendsen barostat does not sample the correct distribution of pressures. 2) Regenerating velocities is fine for a couple of reasons: (a) we are after thermodynamics, not dynamics, and (b) regenerating velocities does not de-equilibrate the system (see the Andersen thermostat, for example). On Jun 8, 2009, at 3:02 PM, Mauricio Carrillo Tripp wrote: Hi, I've been following directions at http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial to calculate hydration free energies of a few small molecules, and I have a couple of questions: 1) Is there a reason for using NVT instead of NPT for the production step? 2) I noticed that between the equilibration and the production runs you generate velocities at each step, instead of using the ones from the previous step. Is there a reason for doing this? wouldn't this 'de-equilibrate' an already equilibrated system? Thank you. -- 0 | Mauricio Carrillo Tripp, PhD / | Department of Molecular Biology, TPC6 0| The Scripps Research Institute \ | 10550 North Torrey Pines Road 0 | La Jolla, California 92037 / | tri...@scripps.edu 0| http://www.scripps.edu/~trippm ** Aut tace aut loquere meliora silentio ** ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
RE: [gmx-users] Re: Free Energy Calculation
A question to the answer below: If the Berendsen barostat does not sample correctly the distribution of pressures then what guarantees that the the volume/density is correct at the end of the equilibration phase when the user is advised to switch to NVT? Peter Nagy The University of Toledo -Original Message- From: gmx-users-boun...@gromacs.org on behalf of David Mobley Sent: Thu 8/6/2009 2:38 PM To: Mauricio Carrillo Tripp Cc: Discussion list for GROMACS users Subject: [gmx-users] Re: Free Energy Calculation Hi, Sorry for the delay answering. These questions are better put on the GROMACS users list. 1) Yes -- the Berendsen barostat does not sample the correct distribution of pressures. 2) Regenerating velocities is fine for a couple of reasons: (a) we are after thermodynamics, not dynamics, and (b) regenerating velocities does not de-equilibrate the system (see the Andersen thermostat, for example). On Jun 8, 2009, at 3:02 PM, Mauricio Carrillo Tripp wrote: Hi, I've been following directions at http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial to calculate hydration free energies of a few small molecules, and I have a couple of questions: 1) Is there a reason for using NVT instead of NPT for the production step? 2) I noticed that between the equilibration and the production runs you generate velocities at each step, instead of using the ones from the previous step. Is there a reason for doing this? wouldn't this 'de-equilibrate' an already equilibrated system? Thank you. -- 0 | Mauricio Carrillo Tripp, PhD / | Department of Molecular Biology, TPC6 0| The Scripps Research Institute \ | 10550 North Torrey Pines Road 0 | La Jolla, California 92037 / | tri...@scripps.edu 0| http://www.scripps.edu/~trippm ** Aut tace aut loquere meliora silentio ** ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Re: Free energy with Gromacs.
You may also want to look at www.alchemistry.org, as it has some general information about how these should be set up. On Wed, Mar 18, 2009 at 8:02 PM, Justin A. Lemkul jalem...@vt.edu wrote: Again, *please keep all Gromacs-related correspondence on the mailing list* The type of information that you have posted is important to share with others who may know how to help. I do not know the answer to your problem. I would suggest contacting the list again with the following information: 1. The relevant contents of your .mdp file(s) 2. What each figure in your image represents (main image vs. inset) -Justin Eudes Fileti wrote: Justin thanks for the reply. Sorry, 16 refers to the value of dg / dl for lambda = 0. I am sending you a link for you see the profile of the curve. http://cromo.ufabc.edu.br/~fileti/web/dgdl-transfer2.jpg This case is for solute in benzene. The peak in lambda=0 is around 3. For sc-power=1 I got a peak much higher than for sc-power=2, and therefore the integration in this case (sc-power=1) could lead me to a much greater error. In both cases (ethanol and benzene) the curves coincide for values of lambda beyond 0.4. bests eef P.S.: I am using the 3.3.3 version. ___ Eudes Eterno Fileti Centro de Ciências Naturais e Humanas Universidade Federal do ABC Rua Santa Adélia, 166 - Bloco B, Sala 1048 09210-170 Santo André - SP Brasil +55.11.4437-8408 skype: eefileti http://cromo.ufabc.edu.br/~fileti/ On Wed, Mar 18, 2009 at 7:13 PM, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: Eudes Fileti wrote: Hello Justin, I am facing a very similar problem to that you experienced and described in (http://www.gromacs.org/pipermail/gmx-users/2008-February/032429.html). I throw this question in the GMX forum and Berk has kindly helped me. But reading the forum I realized that you already could be solved the problem so that maybe I could help more directly. Please keep Gromacs-related discussions on the list. If you had followed the rest of the thread you reference, you will find that my calculations were falling victim to a bug in Gromacs-3.3.1. If you are using a newer version, then my situation is not applicable, since the problem has been fixed (IIRC) as of Gromacs-3.3.3. If you are using version 3.3.1, re-run your simulations with a newer version. I have tried to calculate the free energy of transfer from benzene to ethanol for a polyhydroxylated (24 OH's). This system has 24 hydroxyl groups, and in ethanol, there should be more than 20 solute-solvent hydrogen bonds being erased simultaneously (not to mention the possible intramolecular HB's). The Dg/dlambda plot, for both, benzene and ethanol shows a very high and narrow peak near lambda=0. In the case of ethanol is worse due to the solute-solvent hydrogen bonds. I performed two sets of simulations, one for sc-power=1 and another for sc-power=2, using the following protocol: 1) I made disappear the electrostatic interactions turning off the charges (by 200ps), 2) At the sequence I made disappear the LJ interactions (for more 200ps) 3) Finally I performed a run of 0.5ns. Correct me i this procedure is inappropriate. To start, Berk said me that the use of sc-power=2 never is recommended. Ok! Secondly, then he gave me a good tip that I was not taking into account: Disappear the electrostatic interactions using hard-core instead soft-core. I did this and actually work (only in part). When I used softcore to desappear the electrostatic interations, the value of dg/dlambda for lambda = 0 was ~16. Following the tip of Berk, wiht hardcore I got ~2000! Right, you should only need soft-core for the LJ component. However when I needed to use softcore again, now to remove the LJ interactions, the value returned back to ~16.` I don't understand what you mean. Is the total area under the curve 16? That is ridiculously high. I don't know that I have the expertise to help you much more. I am not entirely familiar with your methodology. What I have done in my own work is run 21 independent simulations at lambda points between 0 and 1 (5 ns each after equilibration), and integrated the resulting curve. None of my dV/dl points ever approached that magnitude, but I can't comment on your specific case, because I don't know what you're doing. -Justin Could you gimme some insigths to solve this problem? Best eef P.S.: You can solve your problem of polyphenolic compound? ___ Eudes Eterno Fileti
[gmx-users] Re: Free energy with Gromacs.
Eudes Fileti wrote: Hello Justin, I am facing a very similar problem to that you experienced and described in (http://www.gromacs.org/pipermail/gmx-users/2008-February/032429.html). I throw this question in the GMX forum and Berk has kindly helped me. But reading the forum I realized that you already could be solved the problem so that maybe I could help more directly. Please keep Gromacs-related discussions on the list. If you had followed the rest of the thread you reference, you will find that my calculations were falling victim to a bug in Gromacs-3.3.1. If you are using a newer version, then my situation is not applicable, since the problem has been fixed (IIRC) as of Gromacs-3.3.3. If you are using version 3.3.1, re-run your simulations with a newer version. I have tried to calculate the free energy of transfer from benzene to ethanol for a polyhydroxylated (24 OH's). This system has 24 hydroxyl groups, and in ethanol, there should be more than 20 solute-solvent hydrogen bonds being erased simultaneously (not to mention the possible intramolecular HB's). The Dg/dlambda plot, for both, benzene and ethanol shows a very high and narrow peak near lambda=0. In the case of ethanol is worse due to the solute-solvent hydrogen bonds. I performed two sets of simulations, one for sc-power=1 and another for sc-power=2, using the following protocol: 1) I made disappear the electrostatic interactions turning off the charges (by 200ps), 2) At the sequence I made disappear the LJ interactions (for more 200ps) 3) Finally I performed a run of 0.5ns. Correct me i this procedure is inappropriate. To start, Berk said me that the use of sc-power=2 never is recommended. Ok! Secondly, then he gave me a good tip that I was not taking into account: Disappear the electrostatic interactions using hard-core instead soft-core. I did this and actually work (only in part). When I used softcore to desappear the electrostatic interations, the value of dg/dlambda for lambda = 0 was ~16. Following the tip of Berk, wiht hardcore I got ~2000! Right, you should only need soft-core for the LJ component. However when I needed to use softcore again, now to remove the LJ interactions, the value returned back to ~16.` I don't understand what you mean. Is the total area under the curve 16? That is ridiculously high. I don't know that I have the expertise to help you much more. I am not entirely familiar with your methodology. What I have done in my own work is run 21 independent simulations at lambda points between 0 and 1 (5 ns each after equilibration), and integrated the resulting curve. None of my dV/dl points ever approached that magnitude, but I can't comment on your specific case, because I don't know what you're doing. -Justin Could you gimme some insigths to solve this problem? Best eef P.S.: You can solve your problem of polyphenolic compound? ___ Eudes Eterno Fileti Centro de Ciências Naturais e Humanas Universidade Federal do ABC Rua Santa Adélia, 166 - Bloco B, Sala 1048 09210-170 Santo André - SP Brasil +55.11.4437-8408 skype: eefileti http://cromo.ufabc.edu.br/~fileti/ -- Justin A. Lemkul Graduate Research Assistant ICTAS Doctoral Scholar 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://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Re: Free energy with Gromacs.
Again, *please keep all Gromacs-related correspondence on the mailing list* The type of information that you have posted is important to share with others who may know how to help. I do not know the answer to your problem. I would suggest contacting the list again with the following information: 1. The relevant contents of your .mdp file(s) 2. What each figure in your image represents (main image vs. inset) -Justin Eudes Fileti wrote: Justin thanks for the reply. Sorry, 16 refers to the value of dg / dl for lambda = 0. I am sending you a link for you see the profile of the curve. http://cromo.ufabc.edu.br/~fileti/web/dgdl-transfer2.jpg This case is for solute in benzene. The peak in lambda=0 is around 3. For sc-power=1 I got a peak much higher than for sc-power=2, and therefore the integration in this case (sc-power=1) could lead me to a much greater error. In both cases (ethanol and benzene) the curves coincide for values of lambda beyond 0.4. bests eef P.S.: I am using the 3.3.3 version. ___ Eudes Eterno Fileti Centro de Ciências Naturais e Humanas Universidade Federal do ABC Rua Santa Adélia, 166 - Bloco B, Sala 1048 09210-170 Santo André - SP Brasil +55.11.4437-8408 skype: eefileti http://cromo.ufabc.edu.br/~fileti/ On Wed, Mar 18, 2009 at 7:13 PM, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: Eudes Fileti wrote: Hello Justin, I am facing a very similar problem to that you experienced and described in (http://www.gromacs.org/pipermail/gmx-users/2008-February/032429.html). I throw this question in the GMX forum and Berk has kindly helped me. But reading the forum I realized that you already could be solved the problem so that maybe I could help more directly. Please keep Gromacs-related discussions on the list. If you had followed the rest of the thread you reference, you will find that my calculations were falling victim to a bug in Gromacs-3.3.1. If you are using a newer version, then my situation is not applicable, since the problem has been fixed (IIRC) as of Gromacs-3.3.3. If you are using version 3.3.1, re-run your simulations with a newer version. I have tried to calculate the free energy of transfer from benzene to ethanol for a polyhydroxylated (24 OH's). This system has 24 hydroxyl groups, and in ethanol, there should be more than 20 solute-solvent hydrogen bonds being erased simultaneously (not to mention the possible intramolecular HB's). The Dg/dlambda plot, for both, benzene and ethanol shows a very high and narrow peak near lambda=0. In the case of ethanol is worse due to the solute-solvent hydrogen bonds. I performed two sets of simulations, one for sc-power=1 and another for sc-power=2, using the following protocol: 1) I made disappear the electrostatic interactions turning off the charges (by 200ps), 2) At the sequence I made disappear the LJ interactions (for more 200ps) 3) Finally I performed a run of 0.5ns. Correct me i this procedure is inappropriate. To start, Berk said me that the use of sc-power=2 never is recommended. Ok! Secondly, then he gave me a good tip that I was not taking into account: Disappear the electrostatic interactions using hard-core instead soft-core. I did this and actually work (only in part). When I used softcore to desappear the electrostatic interations, the value of dg/dlambda for lambda = 0 was ~16. Following the tip of Berk, wiht hardcore I got ~2000! Right, you should only need soft-core for the LJ component. However when I needed to use softcore again, now to remove the LJ interactions, the value returned back to ~16.` I don't understand what you mean. Is the total area under the curve 16? That is ridiculously high. I don't know that I have the expertise to help you much more. I am not entirely familiar with your methodology. What I have done in my own work is run 21 independent simulations at lambda points between 0 and 1 (5 ns each after equilibration), and integrated the resulting curve. None of my dV/dl points ever approached that magnitude, but I can't comment on your specific case, because I don't know what you're doing. -Justin Could you gimme some insigths to solve this problem? Best eef P.S.: You can solve your problem of polyphenolic compound? ___ Eudes Eterno Fileti Centro de Ciências Naturais e Humanas Universidade Federal do ABC Rua Santa Adélia, 166 - Bloco B, Sala 1048 09210-170 Santo André - SP Brasil +55.11.4437-8408 skype:
[gmx-users] Re: free energy with TIP4P bug fixed
Berk, Any chance of getting a fix for 3.3.x versions also? I have several papers which probably are affected by this problem and I will need to repeat the calculations with a fixed version and produce errata. I would prefer to do this with 3.3.x since (a) not all of the data is with TIP4P, and so I don't need to repeat all the calcs, and would like to use a consistent version, and (b) I am not sure that 4.x does not introduce additional bugs that might affect my calcs. Thanks! David On Jan 29, 2009, at 4:21 AM, Berk Hess wrote: Hi, The Coulomb energy difference that Chris Neale observed recently was caused by a bug in the neighborlist assignment with the combination of free energy and tip4p water optimization. This bug would cause a few tip4p-tip4p charge interactions to be missing. I think it has been present in all Gromacs version which have tip4p optimized loops, for sure it was in 3.3. I have fixed this for the upcoming Gromacs 4.0.4 release. I assume this bug also caused the cut-off dependence that David Mobley observed. I have done a lot of free energy calculation with tip4p and never noticed any problems. This was because I always had the perturbed molecule in a separate energy group, which circumvents the problem. So for the moment and for checking if you had the problem with older Gromacs versions, you can simply put the perturbed atoms and tip4p in separate energy groups. Berk What can you do with the new Windows Live? Find out ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
RE: [gmx-users] Re: free energy with TIP4P bug fixed
Hi, You can apply the cvs diff. But even simpler, as I said, is just putting your solute and solvent in separate charge groups. This should circumvent the issue. Berk From: dmob...@gmail.com To: g...@hotmail.com Date: Thu, 29 Jan 2009 06:35:24 -0600 CC: gmx-users@gromacs.org Subject: [gmx-users] Re: free energy with TIP4P bug fixed Berk, Any chance of getting a fix for 3.3.x versions also? I have several papers which probably are affected by this problem and I will need to repeat the calculations with a fixed version and produce errata. I would prefer to do this with 3.3.x since (a) not all of the data is with TIP4P, and so I don't need to repeat all the calcs, and would like to use a consistent version, and (b) I am not sure that 4.x does not introduce additional bugs that might affect my calcs. Thanks!David On Jan 29, 2009, at 4:21 AM, Berk Hess wrote:Hi, The Coulomb energy difference that Chris Neale observed recently was caused by a bug in the neighborlist assignment with the combination of free energy and tip4p water optimization. This bug would cause a few tip4p-tip4p charge interactions to be missing. I think it has been present in all Gromacs version which have tip4p optimized loops, for sure it was in 3.3. I have fixed this for the upcoming Gromacs 4.0.4 release. I assume this bug also caused the cut-off dependence that David Mobley observed. I have done a lot of free energy calculation with tip4p and never noticed any problems. This was because I always had the perturbed molecule in a separate energy group, which circumvents the problem. So for the moment and for checking if you had the problem with older Gromacs versions, you can simply put the perturbed atoms and tip4p in separate energy groups. Berk What can you do with the new Windows Live? Find out _ See all the ways you can stay connected to friends and family http://www.microsoft.com/windows/windowslive/default.aspx___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
RE: [gmx-users] Re: free energy with TIP4P bug fixed
From: g...@hotmail.com To: gmx-users@gromacs.org Subject: RE: [gmx-users] Re: free energy with TIP4P bug fixed Date: Thu, 29 Jan 2009 13:48:20 +0100 Hi, You can apply the cvs diff. But even simpler, as I said, is just putting your solute and solvent in separate charge** groups. ** I meant separate energy groups. This should circumvent the issue. Berk From: dmob...@gmail.com To: g...@hotmail.com Date: Thu, 29 Jan 2009 06:35:24 -0600 CC: gmx-users@gromacs.org Subject: [gmx-users] Re: free energy with TIP4P bug fixed Berk, Any chance of getting a fix for 3.3.x versions also? I have several papers which probably are affected by this problem and I will need to repeat the calculations with a fixed version and produce errata. I would prefer to do this with 3.3.x since (a) not all of the data is with TIP4P, and so I don't need to repeat all the calcs, and would like to use a consistent version, and (b) I am not sure that 4.x does not introduce additional bugs that might affect my calcs. Thanks!David On Jan 29, 2009, at 4:21 AM, Berk Hess wrote:Hi, The Coulomb energy difference that Chris Neale observed recently was caused by a bug in the neighborlist assignment with the combination of free energy and tip4p water optimization. This bug would cause a few tip4p-tip4p charge interactions to be missing. I think it has been present in all Gromacs version which have tip4p optimized loops, for sure it was in 3.3. I have fixed this for the upcoming Gromacs 4.0.4 release. I assume this bug also caused the cut-off dependence that David Mobley observed. I have done a lot of free energy calculation with tip4p and never noticed any problems. This was because I always had the perturbed molecule in a separate energy group, which circumvents the problem. So for the moment and for checking if you had the problem with older Gromacs versions, you can simply put the perturbed atoms and tip4p in separate energy groups. Berk What can you do with the new Windows Live? Find out See all the ways you can stay connected to friends and family _ Express yourself instantly with MSN Messenger! Download today it's FREE! http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Re: Free energy - Hydrogen bond solute-solvent - Ethanol as an example.
Hello all I have tried to reproduce the hydration free energy (TI) of the ethanol from Hess and van der Vegt (JPCB, 110, 17616). The value I have obtained is around 20kJ/mol while the reference value is -20.1kJ/mol (if not the sign ...). If someone can help me find the mistake I would be very grateful. Below are the simulation details. I followed the protocol of the paper, Berk: I used 47 lambda values (because of hydrogen bonds between solute-solvent) (dense near lambda=0 and between 0.46 and 0.72). I am turning off the LJ and Coulomb terms separately. Softcore (alpha = 0.5, power = 1), OPLS-AA Timestep = 2fs, sd = integrator, PME constrained = none trajectories = 40ps (NVT) (I know that is small, but I looking for qualitative results) after 20ps of equilibration (NPT). I have found a value of -4.8kJ/mol for the DeltaG(vacuum) (relative to mutate the ethanol to dummy in vacuum). My dv/dl curve (for DeltaG(water) is below. Should I expect this form? I think there are many high positive values. The numerical integration is 15.9kJ. So, DG(hyd) = DG(wat) - DG(vac) = 15.9 - (-4.8) = 20.7kJ I believe that this protocol is OK (but I want to confirm that). @title dG/d\8l\4 @xaxis label Time (ps) @yaxis label dG/d\8l\4 (kJ mol\S-1\N [\8l\4]\S-1\N) @TYPE xy 0.000 4.431073e+02 0.005 3.660083e+02 0.010 2.987275e+02 0.015 2.023206e+02 0.000 4.431073e+02 0.005 3.660083e+02 0.010 2.987275e+02 0.015 2.023206e+02 0.020 1.905198e+02 0.030 1.086573e+02 0.040 4.217960e+01 0.050 4.101856e+01 0.060 2.314358e+01 0.070 2.223774e+01 0.080 2.669019e+01 0.090 1.696131e+01 0.100 6.089735e+00 0.110 9.569030e+00 0.120 1.505562e+01 0.130 2.785974e+00 0.140 2.440906e+00 0.150 3.050463e+00 0.160 8.578429e-01 0.200 3.577259e+00 0.240 -3.532969e+00 0.280 5.842623e+00 0.320 7.912565e+00 0.360 9.322502e+00 0.400 -2.930754e+00 0.440 2.797944e+00 0.460 5.113900e+00 0.480 -3.919665e+00 0.500 -9.000282e+00 0.520 -7.636168e+00 0.540 -7.574299e+00 0.560 -2.144124e+01 0.580 -6.618514e+00 0.600 -2.755925e+01 0.620 -1.491331e+01 0.640 -2.541398e+01 0.660 -5.213154e+01 0.680 -2.364612e+01 0.700 -2.535817e+01 0.720 -7.767502e+00 0.760 -1.096362e+01 0.800 1.136031e+01 0.840 2.883933e+01 0.880 4.465743e+01 0.920 5.461566e+01 0.960 6.882000e+01 1.000 7.673736e+01 -- ___ Eudes Eterno Fileti Centro de Ciência Naturais e Humanas Universidade Federal do ABC Rua Catequese, 242 - 3º Andar 09090-400 Santo André - SP Brasil Tel: +55 11 4437-1600 ramal 408 skype: eefileti http://cromo.ufabc.edu.br/~fileti/ ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Re: Free energy
Dear Eudes, I don't have time to provide personalized help with this at present. You'll be better off writing the gromacs forum. You may also want to refer to the general information on performing free energy calculations at http://www.alchemistry.org and previous posts on the gromacs users list. I have personally done ethanol hydration free energies myself with no problems, so if you are getting wrong values it is probably because you have done something wrong -- most likely a question I have already answered on the users list or at alchemistry.org. Best wishes, David On Tue, Apr 8, 2008 at 11:28 AM, Eudes Fileti [EMAIL PROTECTED] wrote: Dear Mobley I'm sorry for writing you (instead to gmx-forum). I have read your tutorial on the free energy for methane and I reproduced the results that you presented there. However, now I am trying a slightly more complicated system (ethanol) and I can not performed the calculation. The results that I am getting, (using the same protocol that you used at JPCB 111,2242) are very larger than the typical values (on the order of 10E+05). Moreover the simulations only conclude for values of init_lamba lower than 0.60, for higher values, the simulation chashes with segmentation fault) I send you my .top and .mdp files. If you can help me I will very grateful Thanks eef -- ___ Eudes Eterno Fileti Centro de Ciência Naturais e Humanas Universidade Federal do ABC Rua Catequese, 242 - 3º Andar 09090-400 Santo André - SP Brasil Tel: +55 11 4437-1600 ramal 408 skype: eefileti http://cromo.ufabc.edu.br/~fileti/ ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Re: free energy calculation
Hi, I am traveling and don't have time to answer this now, nor am I responsible for the tutorial you mention. I also only usually answer questions like this on the gromacs users list, not in my personal e-mail. I suggest you direct your e-mail to the Gromacs users list (which I am ccing on this) and get help from someone else on there, or i will be able to attend to it (in that forum) in a few days. Thanks, David Mobley On Feb 11, 2008 5:14 AM, [EMAIL PROTECTED] wrote: Dear david, i 'm a beginner in molecular dynamics. I was doing the tutorial of md group, hydration free energy of toluene : http://md.chem.rug.nl/education/Free-Energy_Course/2.hydration-fe.html i have some problems regarding the calculation of toluene in water. the thing is that in the examples files equ-10.00.mdp and data-10.00.mdp no sc-power were mentioned and doing the calculation the error files mentioned about sc-power must be 0 so i added this line sc-power=1 to those files. Analyzing the dgdl.xvg files the estimate error(err. est.) is about 5.79 for each lambda value studied while the mutation of toluene to dummy in vacuo gave really reasonable rms of about 10E-5. Combining the integration value of the plot in water and vacuo did not give me the right solvation free energy value (compare to experimental one). i assumed that this is the free calculation in water which is problematic. I tried to play with sc-power and sc-alpha but still even if those parameters seem to influence a lot on the energy no real improvement are visible. Would you have any idea what could be the problem? geraldine Helsinki university ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] RE: free energy calculation
Dear all, Can anybody tell me how to calculate the free energy difference after running simulation with lamda(0, 0.1,...1)? I am following the tutorial on wiki, but i have no idea about the details to do the data analysis. any software can do it? and any reference for the instruction? I am new here. thanks a lot Qiang Try to read this: http://en.wikipedia.org/wiki/Numerical_integration Vojtech Spiwok winmail.dat___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Re: free energy tutorial
David Mobley wrote: Bharat, Matt, Tsjerk, and all, I just put a bit of a tutorial on free energy calculations up at http://www.dillgroup.ucsf.edu/~jchodera/group/wiki/index.php/Free_Energy:_Tutorial. This will probably be a work in progress; it is not as step-by-step as Tsjerk's (although possibly I could make it that way later), but it IS current for GROMACS 3.3/3.3.1. Feel free to contact me directly if there are any comments/questions/problems. I should mention at this point the tutorial just covers disappearing neutral methane in water; I'll probably expand it later. I do think this should probably be a standard test case for people starting off on free energy calculations, as Michael Shirts and I have both calculated values for this in GROMACS with OPLS-AA (and our results agree quite well, always a plus). Thanks David, I've now put up a list of on-line tutorials under Documentation/Tutorial on the website. Please let me know if there are other resources that could be added. -- David. David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group, Dept. of Cell and Molecular Biology, Uppsala University. Husargatan 3, Box 596, 75124 Uppsala, Sweden phone: 46 18 471 4205 fax: 46 18 511 755 [EMAIL PROTECTED] [EMAIL PROTECTED] http://folding.bmc.uu.se ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Re: free energy tutorial
Bharat, Matt, Tsjerk, and all, I just put a bit of a tutorial on free energy calculations up at http://www.dillgroup.ucsf.edu/~jchodera/group/wiki/index.php/Free_Energy:_Tutorial. This will probably be a work in progress; it is not as step-by-step as Tsjerk's (although possibly I could make it that way later), but it IS current for GROMACS 3.3/3.3.1. Feel free to contact me directly if there are any comments/questions/problems. I should mention at this point the tutorial just covers disappearing neutral methane in water; I'll probably expand it later. I do think this should probably be a standard test case for people starting off on free energy calculations, as Michael Shirts and I have both calculated values for this in GROMACS with OPLS-AA (and our results agree quite well, always a plus). Thanks, David Mobley UCSF On 6/4/06, Tsjerk Wassenaar [EMAIL PROTECTED] wrote: Hi Bharat and Matt, The tutorial was set up by a colleague of mine, who has left our group a while ago. Given your remarks, I think it's better to replace the current version with a new one. I hope I can find some time for the revision and would welcome any suggestions for changes, additions (topics to cover), which can be sent to me off the list. Consider the tutorial in its current form deprecated (which I will add as a note in the tutorial soon). Best regards, Tsjerk On 6/2/06, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: Hi Bharat, This is the same tutorial that led me astray when I first began trying free energy calculations in Gromacs. I am not sure how close to the experimental value the author intended to get with this hydration free energy tutorial. If you are using a recent version of Gromacs, be aware that the soft core method has changed. The tutorial says sc-alpha should be 1.51, but that is more appropriate for the older method corresponding to sc-power=2 rather than the sc-power=1 that is now recommended. This tutorial also has you simultaneously mutate partial charges and atom types. I have recently been informed, and confirmed via observation, that performing the partial charge mutation with soft core potentials on can lead to misleading and very noisy results. Look in the recent mailing list archives for suggestions and papers about splitting the transformation into two series of calculations, one with partial charge mutation and no soft core potentials and the other with soft core potentials and no partial charge mutation. Even if you don't want to do this for the simple tutorial you'll probably want to do it for real systems in the future. Finally, the tutorial gives the dubious advice to use the final structure from one window as the starting structure for the next window. If you instead run each window's simulations completely independently, you can simultaneously run calculations on several machines and don't have to worry about the output of previous windows somehow biasing or contaminating the current window. Again, I don't know if it is necessary to operate this way to get good results on this simple tutorial, but it does become more important if you want to run novel, non-trivial systems in the future. Matt Ernst Washington State University Dear GMX users, i am trying 'relative hydration free energy of p-cresol wrt toluene' tutorial from http://md.chem.rug.nl/education/Free-Energy_Course/3.mutation-fe.html . i am using all the mdp files given on the site, the confout.gro from toluene/eq directory of the previous tutorial as starting structure for this simulation. and whatever relevant changes has to be made in top and gro have been made. then i performed the simulation for lamda values from zero to one with an interval of 0.4, so total of 26 simulations. but i am not getting the mentioned experimental value of ddG-hyd. i got the following values in rmsd.xvg IN WATERIN VACUO == 0-61.874991 2.591030-98.521845 8.80653 .04 -57.180709 1.73975.04 -67.562344 0.237903 .08 -58.786951 0.266783 .08 -64.780817 1.23656 .12 -60.764221 0.300319 .12 -63.280561 1.14712 .16 -62.7007740.313339 .16 -60.123999 0.926409 .20 -65.107642 1.26388.20 -58.584394 0.103875 .24 -63.749578 0.384702 .24 -56.668703 1.02234 .28 -64.267572 1.22184.28 -54.837984 0.0807904 .32 -61.524424 0.402878 .32 -53.550764 0.0760657 .36 -55.897929 0.314078 .36 -51.510805 0.0622261 .40 -49.084890 1.05137.40 -49.736461 0.0627338 .44 -38.3398011.39297.44 -48.470839 0.0718284 .48 -29.613058 1.38755.48 -46.842750 0.0740531 .52 -24.891702 1.04125.52 -45.007666 0.0766915 .56 -25.435188 1.70549.56 -43.9396820.101024 .60 -25.180946 1.82089.60 -41.844138 0.111544 .64 -28.395296 0.849711 .64 -41.250909 0.13644 .68
Re: [gmx-users] Re: free energy tutorial
Hi Bharat and Matt, The tutorial was set up by a colleague of mine, who has left our group a while ago. Given your remarks, I think it's better to replace the current version with a new one. I hope I can find some time for the revision and would welcome any suggestions for changes, additions (topics to cover), which can be sent to me off the list. Consider the tutorial in its current form deprecated (which I will add as a note in the tutorial soon). Best regards, TsjerkOn 6/2/06, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote:Hi Bharat,This is the same tutorial that led me astray when I first began trying free energy calculations in Gromacs. I am not sure how close to the experimental value the authorintended to get with this hydration free energy tutorial. If you are using a recentversion of Gromacs, be aware that the soft core method has changed. The tutorial says sc-alpha should be 1.51, but that is more appropriate for the older method correspondingto sc-power=2 rather than the sc-power=1 that is now recommended.This tutorial also has you simultaneously mutate partial charges and atom types. I have recently been informed, and confirmed via observation, that performing the partialcharge mutation with soft core potentials on can lead to misleading and very noisyresults. Look in the recent mailing list archives for suggestions and papers about splitting the transformation into two series of calculations, one with partial chargemutation and no soft core potentials and the other with soft core potentials and nopartial charge mutation. Even if you don't want to do this for the simple tutorial you'll probably want to do it for real systems in the future.Finally, the tutorial gives the dubious advice to use the final structure from onewindow as the starting structure for the next window. If you instead run each window's simulations completely independently, you can simultaneously run calculations on severalmachines and don't have to worry about the output of previous windows somehow biasing orcontaminating the current window. Again, I don't know if it is necessary to operate this way to get good results on this simple tutorial, but it does become more important ifyou want to run novel, non-trivial systems in the future.Matt ErnstWashington State University Dear GMX users, i am trying 'relative hydration free energy of p-cresol wrt toluene' tutorial from http://md.chem.rug.nl/education/Free-Energy_Course/3.mutation-fe.html . i am using all the mdp files given on the site, the confout.gro from toluene/eq directory of the previous tutorial as starting structure for this simulation. and whatever relevant changes has to be made in top and gro have been made. then i performed the simulation for lamda values from zero to one with an interval of 0.4, so total of 26 simulations. but i am not getting the mentioned experimental value of ddG-hyd. i got the following values in rmsd.xvg IN WATERIN VACUO == 0-61.874991 2.591030-98.521845 8.80653 .04-57.180709 1.73975.04-67.562344 0.237903 .08-58.786951 0.266783 .08-64.780817 1.23656 .12-60.764221 0.300319 .12-63.280561 1.14712 .16-62.700774 0.313339 .16-60.123999 0.926409 .20-65.107642 1.26388.20-58.584394 0.103875 .24-63.749578 0.384702 .24-56.668703 1.02234 .28-64.267572 1.22184.28-54.837984 0.0807904 .32-61.524424 0.402878 .32-53.550764 0.0760657 .36-55.897929 0.314078 .36-51.510805 0.0622261 .40-49.084890 1.05137.40-49.736461 0.0627338 .44-38.339801 1.39297.44-48.470839 0.0718284 .48-29.613058 1.38755.48-46.842750 0.0740531 .52-24.891702 1.04125.52-45.007666 0.0766915 .56-25.435188 1.70549.56-43.939682 0.101024 .60-25.180946 1.82089.60-41.844138 0.111544 .64-28.395296 0.849711 .64-41.250909 0.13644 .68-29.672103 2.26016.68-41.182423 0.139366 .72-34.441109 2.00941.72-41.747919 0.127424 .76-39.148106 2.01051.76-44.426396 0.143008 .80-49.360365 2.34239.80-44.381975 0.135462 .84-57.031557 2.79298.84-45.564523 0.123685 .88-64.059494 2.79089.88-46.074452 0.104784 .92-78.668957 1.97957.92-46.678656 0.0993529 .96-85.671014 1.92346.96-46.871521 0.0923201 1.0-52.956186 1.261361.0-47.183714 0.0934921 upon integration in xmgr, following are the dG values: dG-wat = -50.6556 kJ.mol-1 dG-vac = -51.1109 kJ.mol-1 ddG-hyd = 0.455 kJ.mol-1 the experimental value mentioned on the site is -22.1 kJ.mol-1 Please help... bharat___gmx-users mailing list gmx-users@gromacs.orghttp://www.gromacs.org/mailman/listinfo/gmx-usersPlease don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED].Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- Tsjerk A. Wassenaar, M.Sc.Groningen Biomolecular Sciences and Biotechnology Institute (GBB)Dept. of Biophysical ChemistryUniversity of Groningen Nijenborgh 49747AG Groningen, The Netherlands+31 50 363 4336 ___ gmx-users mailing listgmx-users@gromacs.org
[gmx-users] Re: free energy tutorial
Hi Bharat, This is the same tutorial that led me astray when I first began trying free energy calculations in Gromacs. I am not sure how close to the experimental value the author intended to get with this hydration free energy tutorial. If you are using a recent version of Gromacs, be aware that the soft core method has changed. The tutorial says sc-alpha should be 1.51, but that is more appropriate for the older method corresponding to sc-power=2 rather than the sc-power=1 that is now recommended. This tutorial also has you simultaneously mutate partial charges and atom types. I have recently been informed, and confirmed via observation, that performing the partial charge mutation with soft core potentials on can lead to misleading and very noisy results. Look in the recent mailing list archives for suggestions and papers about splitting the transformation into two series of calculations, one with partial charge mutation and no soft core potentials and the other with soft core potentials and no partial charge mutation. Even if you don't want to do this for the simple tutorial you'll probably want to do it for real systems in the future. Finally, the tutorial gives the dubious advice to use the final structure from one window as the starting structure for the next window. If you instead run each window's simulations completely independently, you can simultaneously run calculations on several machines and don't have to worry about the output of previous windows somehow biasing or contaminating the current window. Again, I don't know if it is necessary to operate this way to get good results on this simple tutorial, but it does become more important if you want to run novel, non-trivial systems in the future. Matt Ernst Washington State University Dear GMX users, i am trying 'relative hydration free energy of p-cresol wrt toluene' tutorial from http://md.chem.rug.nl/education/Free-Energy_Course/3.mutation-fe.html. i am using all the mdp files given on the site, the confout.gro from toluene/eq directory of the previous tutorial as starting structure for this simulation. and whatever relevant changes has to be made in top and gro have been made. then i performed the simulation for lamda values from zero to one with an interval of 0.4, so total of 26 simulations. but i am not getting the mentioned experimental value of ddG-hyd. i got the following values in rmsd.xvg IN WATERIN VACUO == 0-61.874991 2.591030-98.521845 8.80653 .04 -57.180709 1.73975.04 -67.562344 0.237903 .08 -58.786951 0.266783 .08 -64.780817 1.23656 .12 -60.764221 0.300319 .12 -63.280561 1.14712 .16 -62.700774 0.313339 .16 -60.123999 0.926409 .20 -65.107642 1.26388.20 -58.584394 0.103875 .24 -63.749578 0.384702 .24 -56.668703 1.02234 .28 -64.267572 1.22184.28 -54.837984 0.0807904 .32 -61.524424 0.402878 .32 -53.550764 0.0760657 .36 -55.897929 0.314078 .36 -51.510805 0.0622261 .40 -49.084890 1.05137.40 -49.736461 0.0627338 .44 -38.339801 1.39297.44 -48.470839 0.0718284 .48 -29.613058 1.38755.48 -46.842750 0.0740531 .52 -24.891702 1.04125.52 -45.007666 0.0766915 .56 -25.435188 1.70549.56 -43.939682 0.101024 .60 -25.180946 1.82089.60 -41.844138 0.111544 .64 -28.395296 0.849711 .64 -41.250909 0.13644 .68 -29.672103 2.26016.68 -41.182423 0.139366 .72 -34.441109 2.00941.72 -41.747919 0.127424 .76 -39.148106 2.01051.76 -44.426396 0.143008 .80 -49.360365 2.34239.80 -44.381975 0.135462 .84 -57.031557 2.79298.84 -45.564523 0.123685 .88 -64.059494 2.79089.88 -46.074452 0.104784 .92 -78.668957 1.97957.92 -46.678656 0.0993529 .96 -85.671014 1.92346.96 -46.871521 0.0923201 1.0 -52.956186 1.261361.0 -47.183714 0.0934921 upon integration in xmgr, following are the dG values: dG-wat = -50.6556 kJ.mol-1 dG-vac = -51.1109 kJ.mol-1 ddG-hyd = 0.455 kJ.mol-1 the experimental value mentioned on the site is -22.1 kJ.mol-1 Please help... bharat ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
