Re: [gmx-users] Dubious results with NPT
Hello, This time I tried sesiisotropic scaling with Berendsen Barostat. The compressibility was set to zero in Z-direction allowing scaling only in X and Y. The NPT simulation was run for 10ns. The result did not look good and the average pressure was ~36 bar (reference pressure is 1 bar). Even more confusing was the graphs of various properties against time. Pressure vs time shows a sudden decrease in fluctuations at about 5ns and thereafter . The pressure fluctuations lowered to around ~170bar from ~650 bar. A dramatic change was observed in many properties at around this time. For eg, the short range LJ interaction was found to decrease sharply at ~5ns. Also, total energy, Kinetic, Potential energies and the temperature showed a rapid drop in the magnitude of their fluctuations from the same time. Nothing notable was observed in the system trajectory at around this time. I am totally confused, and I feel that I made some very bad choice again the MDP file given below. Please let me know what you think. Suggest something on this issue which I should read about. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.001 nsteps = 1000 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 25000 nstvout = 25000 nstfout = 25000 nstlog = 100 nstenergy = 100 nstxtcout = 100 xtc_precision = 1000 xtc_grps = System energygrps = System ; 7.3.9 Neighbor Searching nstlist= 10 ns_type = grid pbc = xyz rlist = 0.9 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.9 ; 7.3.11 VdW vdwtype = cut-off rvdw = 1.4 DispCorr = EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol= 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t = 0.50.5 ref_t = 298298 ; 7.3.15 Pressure Coupling pcoupl = berendsen pcoupltype = semiisotropic tau_p = 2.0 compressibility= 4.5e-5 0.0 ref_p = 1.01.0 ; 7.3.17 Velocity Generation gen_vel = no gen_temp= 298 gen_seed = -1 ; 7.3.18 Bonds constraints = none constraint_algorithm= LINCS continuation= no lincs_order = 4 lincs_iter = 1 lincs_warnangle = 30 Regards, Sujith. On Wed, Mar 12, 2014 at 5:28 PM, Justin Lemkul jalem...@vt.edu wrote: On 3/12/14, 7:51 AM, sujithkakkat . wrote: Hello, After a while I got back to the problem posted here. This issue was the large value for average pressure(~25 bar against the reference pressure of 1 bar) in NPT simulations with parrinello rahman barostat. The system studied is cyclohexane-water system with an interface. The forcefield is Gromos96. Gromacs 4.6.5 version is used. I recollect a post from Michael Shirts which says that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. Is this what is happening in my case? Based on your suggestions , I tried the following to analyze/solve the problem, and till now the results have not improved a bit. (i) The system was simplified. Each phases was simulated separately. The results were good, with proper values of average pressure after NPT simulation for both independent water and cyclohexane systems. However, when doing the same with similar parameters on the water-cyclohexane combined system with parrinello-rahman barostat , still gives very high average pressure (even after 10ns ) (ii) Changing the tau_p values did not help either. (I tried 1ps, 2ps, and 4ps) My question is whether, I can try the semiisotropic scaling? The manual says that it is useful in case of systems with interface. I am not sure, since none of the replies I got suggested that option, which gives the impression that it is not the best solution. But , to my limited knowledge that is the only option left to be tried to solve the issue. Worth a shot. Also try changing the barostat to Berendsen to see if that changes anything. There are a lot of interrelated issues here, possibly a problem with the barostat, the tau_p value, the type of coupling, the combination of the integrator and barostat...you see how complicated it gets, and that's only some of the possible issues. -Justin The mdp file used is, ; 7.3.3 Run Control integrator = md
Re: [gmx-users] Dubious results with NPT
The simplest explanation would be that you've appended to a previous 5ns trajectory, not run a new trajectory. Check the .log file and the length of time you expected this job to run (wall and simulation). Mark On Fri, Mar 14, 2014 at 9:02 AM, sujithkakkat . sujithk...@gmail.comwrote: Hello, This time I tried sesiisotropic scaling with Berendsen Barostat. The compressibility was set to zero in Z-direction allowing scaling only in X and Y. The NPT simulation was run for 10ns. The result did not look good and the average pressure was ~36 bar (reference pressure is 1 bar). Even more confusing was the graphs of various properties against time. Pressure vs time shows a sudden decrease in fluctuations at about 5ns and thereafter . The pressure fluctuations lowered to around ~170bar from ~650 bar. A dramatic change was observed in many properties at around this time. For eg, the short range LJ interaction was found to decrease sharply at ~5ns. Also, total energy, Kinetic, Potential energies and the temperature showed a rapid drop in the magnitude of their fluctuations from the same time. Nothing notable was observed in the system trajectory at around this time. I am totally confused, and I feel that I made some very bad choice again the MDP file given below. Please let me know what you think. Suggest something on this issue which I should read about. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.001 nsteps = 1000 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 25000 nstvout = 25000 nstfout = 25000 nstlog = 100 nstenergy = 100 nstxtcout = 100 xtc_precision = 1000 xtc_grps = System energygrps = System ; 7.3.9 Neighbor Searching nstlist= 10 ns_type = grid pbc = xyz rlist = 0.9 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.9 ; 7.3.11 VdW vdwtype = cut-off rvdw = 1.4 DispCorr = EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol= 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t = 0.50.5 ref_t = 298298 ; 7.3.15 Pressure Coupling pcoupl = berendsen pcoupltype = semiisotropic tau_p = 2.0 compressibility= 4.5e-5 0.0 ref_p = 1.01.0 ; 7.3.17 Velocity Generation gen_vel = no gen_temp= 298 gen_seed = -1 ; 7.3.18 Bonds constraints = none constraint_algorithm= LINCS continuation= no lincs_order = 4 lincs_iter = 1 lincs_warnangle = 30 Regards, Sujith. On Wed, Mar 12, 2014 at 5:28 PM, Justin Lemkul jalem...@vt.edu wrote: On 3/12/14, 7:51 AM, sujithkakkat . wrote: Hello, After a while I got back to the problem posted here. This issue was the large value for average pressure(~25 bar against the reference pressure of 1 bar) in NPT simulations with parrinello rahman barostat. The system studied is cyclohexane-water system with an interface. The forcefield is Gromos96. Gromacs 4.6.5 version is used. I recollect a post from Michael Shirts which says that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. Is this what is happening in my case? Based on your suggestions , I tried the following to analyze/solve the problem, and till now the results have not improved a bit. (i) The system was simplified. Each phases was simulated separately. The results were good, with proper values of average pressure after NPT simulation for both independent water and cyclohexane systems. However, when doing the same with similar parameters on the water-cyclohexane combined system with parrinello-rahman barostat , still gives very high average pressure (even after 10ns ) (ii) Changing the tau_p values did not help either. (I tried 1ps, 2ps, and 4ps) My question is whether, I can try the semiisotropic scaling? The manual says that it is useful in case of systems with interface. I am not sure, since none of the replies I got suggested that option, which gives the impression that it is not the best solution. But , to my limited knowledge that is the only option left to be tried to solve the issue. Worth a shot. Also try changing the
Re: [gmx-users] Dubious results with NPT
Hello Mark, I guess you were asking whether I ran the simulation as a continuation to a previous 5ns run? But the results I got are from a continuous 10ns simulation. Sujith. On Fri, Mar 14, 2014 at 4:28 PM, Mark Abraham mark.j.abra...@gmail.comwrote: The simplest explanation would be that you've appended to a previous 5ns trajectory, not run a new trajectory. Check the .log file and the length of time you expected this job to run (wall and simulation). Mark On Fri, Mar 14, 2014 at 9:02 AM, sujithkakkat . sujithk...@gmail.com wrote: Hello, This time I tried sesiisotropic scaling with Berendsen Barostat. The compressibility was set to zero in Z-direction allowing scaling only in X and Y. The NPT simulation was run for 10ns. The result did not look good and the average pressure was ~36 bar (reference pressure is 1 bar). Even more confusing was the graphs of various properties against time. Pressure vs time shows a sudden decrease in fluctuations at about 5ns and thereafter . The pressure fluctuations lowered to around ~170bar from ~650 bar. A dramatic change was observed in many properties at around this time. For eg, the short range LJ interaction was found to decrease sharply at ~5ns. Also, total energy, Kinetic, Potential energies and the temperature showed a rapid drop in the magnitude of their fluctuations from the same time. Nothing notable was observed in the system trajectory at around this time. I am totally confused, and I feel that I made some very bad choice again the MDP file given below. Please let me know what you think. Suggest something on this issue which I should read about. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.001 nsteps = 1000 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 25000 nstvout = 25000 nstfout = 25000 nstlog = 100 nstenergy = 100 nstxtcout = 100 xtc_precision = 1000 xtc_grps = System energygrps = System ; 7.3.9 Neighbor Searching nstlist= 10 ns_type = grid pbc = xyz rlist = 0.9 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.9 ; 7.3.11 VdW vdwtype = cut-off rvdw = 1.4 DispCorr = EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol= 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t = 0.50.5 ref_t = 298298 ; 7.3.15 Pressure Coupling pcoupl = berendsen pcoupltype = semiisotropic tau_p = 2.0 compressibility= 4.5e-5 0.0 ref_p = 1.01.0 ; 7.3.17 Velocity Generation gen_vel = no gen_temp= 298 gen_seed = -1 ; 7.3.18 Bonds constraints = none constraint_algorithm= LINCS continuation= no lincs_order = 4 lincs_iter = 1 lincs_warnangle = 30 Regards, Sujith. On Wed, Mar 12, 2014 at 5:28 PM, Justin Lemkul jalem...@vt.edu wrote: On 3/12/14, 7:51 AM, sujithkakkat . wrote: Hello, After a while I got back to the problem posted here. This issue was the large value for average pressure(~25 bar against the reference pressure of 1 bar) in NPT simulations with parrinello rahman barostat. The system studied is cyclohexane-water system with an interface. The forcefield is Gromos96. Gromacs 4.6.5 version is used. I recollect a post from Michael Shirts which says that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. Is this what is happening in my case? Based on your suggestions , I tried the following to analyze/solve the problem, and till now the results have not improved a bit. (i) The system was simplified. Each phases was simulated separately. The results were good, with proper values of average pressure after NPT simulation for both independent water and cyclohexane systems. However, when doing the same with similar parameters on the water-cyclohexane combined system with parrinello-rahman barostat , still gives very high average pressure (even after 10ns ) (ii) Changing the tau_p values did not help either. (I tried 1ps, 2ps, and 4ps)
Re: [gmx-users] Dubious results with NPT
Hello, After a while I got back to the problem posted here. This issue was the large value for average pressure(~25 bar against the reference pressure of 1 bar) in NPT simulations with parrinello rahman barostat. The system studied is cyclohexane-water system with an interface. The forcefield is Gromos96. Gromacs 4.6.5 version is used. I recollect a post from Michael Shirts which says that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. Is this what is happening in my case? Based on your suggestions , I tried the following to analyze/solve the problem, and till now the results have not improved a bit. (i) The system was simplified. Each phases was simulated separately. The results were good, with proper values of average pressure after NPT simulation for both independent water and cyclohexane systems. However, when doing the same with similar parameters on the water-cyclohexane combined system with parrinello-rahman barostat , still gives very high average pressure (even after 10ns ) (ii) Changing the tau_p values did not help either. (I tried 1ps, 2ps, and 4ps) My question is whether, I can try the semiisotropic scaling? The manual says that it is useful in case of systems with interface. I am not sure, since none of the replies I got suggested that option, which gives the impression that it is not the best solution. But , to my limited knowledge that is the only option left to be tried to solve the issue. The mdp file used is, ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.001 nsteps = 1000 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 25000 nstvout = 25000 nstfout = 25000 nstlog = 100 nstenergy = 100 nstxtcout = 100 xtc_precision = 1000 xtc_grps = System energygrps = System ; 7.3.9 Neighbor Searching nstlist= 10 ns_type = grid pbc = xyz rlist = 0.9 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.9 ; 7.3.11 VdW vdwtype = cut-off rvdw = 1.4 DispCorr = EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol= 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t = 0.50.5 ref_t = 298298 ; 7.3.15 Pressure Coupling pcoupl = parrinello-rahman pcoupltype = isotropic tau_p = 2.0 compressibility= 4.5e-5 ref_p = 1.0 ; 7.3.17 Velocity Generation gen_vel = no gen_temp= 298 gen_seed = -1 ; 7.3.18 Bonds constraints = none constraint_algorithm= LINCS continuation= no lincs_order = 4 lincs_iter = 1 lincs_warnangle = 30 Please comment. Sujith. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Dubious results with NPT
On 3/12/14, 7:51 AM, sujithkakkat . wrote: Hello, After a while I got back to the problem posted here. This issue was the large value for average pressure(~25 bar against the reference pressure of 1 bar) in NPT simulations with parrinello rahman barostat. The system studied is cyclohexane-water system with an interface. The forcefield is Gromos96. Gromacs 4.6.5 version is used. I recollect a post from Michael Shirts which says that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. Is this what is happening in my case? Based on your suggestions , I tried the following to analyze/solve the problem, and till now the results have not improved a bit. (i) The system was simplified. Each phases was simulated separately. The results were good, with proper values of average pressure after NPT simulation for both independent water and cyclohexane systems. However, when doing the same with similar parameters on the water-cyclohexane combined system with parrinello-rahman barostat , still gives very high average pressure (even after 10ns ) (ii) Changing the tau_p values did not help either. (I tried 1ps, 2ps, and 4ps) My question is whether, I can try the semiisotropic scaling? The manual says that it is useful in case of systems with interface. I am not sure, since none of the replies I got suggested that option, which gives the impression that it is not the best solution. But , to my limited knowledge that is the only option left to be tried to solve the issue. Worth a shot. Also try changing the barostat to Berendsen to see if that changes anything. There are a lot of interrelated issues here, possibly a problem with the barostat, the tau_p value, the type of coupling, the combination of the integrator and barostat...you see how complicated it gets, and that's only some of the possible issues. -Justin The mdp file used is, ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.001 nsteps = 1000 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 25000 nstvout = 25000 nstfout = 25000 nstlog = 100 nstenergy = 100 nstxtcout = 100 xtc_precision = 1000 xtc_grps = System energygrps = System ; 7.3.9 Neighbor Searching nstlist= 10 ns_type = grid pbc = xyz rlist = 0.9 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.9 ; 7.3.11 VdW vdwtype = cut-off rvdw = 1.4 DispCorr = EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol= 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t = 0.50.5 ref_t = 298298 ; 7.3.15 Pressure Coupling pcoupl = parrinello-rahman pcoupltype = isotropic tau_p = 2.0 compressibility= 4.5e-5 ref_p = 1.0 ; 7.3.17 Velocity Generation gen_vel = no gen_temp= 298 gen_seed = -1 ; 7.3.18 Bonds constraints = none constraint_algorithm= LINCS continuation= no lincs_order = 4 lincs_iter = 1 lincs_warnangle = 30 Please comment. Sujith. -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 601 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Dubious results with NPT
your average pressure is the pressure you should report in the publication. if you got 5 bars instead of 1 bar, you should write I simulated the system at 5 bars Dr. Vitaly V. Chaban On Tue, Feb 25, 2014 at 6:01 AM, sujithkakkat . sujithk...@gmail.com wrote: On Mon, Feb 24, 2014 at 5:37 AM, sujithkakkat . sujithk...@gmail.com wrote: Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. No, this average pressure is not fine. Why isn't it fine? I admit it is not very good. I guess due to huge fluctuations the average pressure is tend to vary a bit to either side of the reference value, unless I continue the equilibration for a very long time. Also the average was calculated with g_energy, which I guess considered the entire simulation, where the first few time steps, when the system is away from equilibrium, would have poor pressure values with respect to reference , which would reflect in the average pressure calculated. Also using g_analyze to calculate average pressure between different times during the simulation did not give the same value. So I am not sure if I should aim at an average pressure which is equal to the reference pressure. Thanks, Sujith. On Mon, Feb 24, 2014 at 4:26 PM, Justin Lemkul jalem...@vt.edu wrote: On 2/23/14, 11:37 PM, sujithkakkat . wrote: Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. I guess the sigma for the case of water was low and therefore the small cut-off of 0.8nm did not matter. However the case of cyclohexane alone remains to be tried. I guess Dr Vitaly was saying about using a switch/shift function. I will try the simulation with the new settings and see. The Gromos force fields were parametrized with a plain cutoff for van der Waals. They do not require switching or shifting. It's true that plain cutoffs are a bit harsh, but switching and shifting have their own issues, and one should always stick to the parametrization methods unless there is demonstrable evidence that the new method is better. I have never known Gromos96 force fields to need anything other than plain cutoff for van der Waals to function as expected. -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 601 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Dubious results with NPT
No full solutions to this problem -- I'll write a few notes. * With large systems run for relatively small amounts of time, the ensemble could be statistically indistinguishble from the true distribution even if the averages don't line up correctly. See: http://pubs.acs.org/doi/abs/10.1021/ct300688p For ways to examine the full distribution of volumes. * Parrinello-Rahman isn't quite correct in GROMACS because virial and kinetic energy are offset by 1/2 a step, though getting something more like 1.2 atm when running at 1 atm is more common. * Note that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. * If you have a property that changes when the pressure changes from 1 to 5 atm, odds are that molecular simulation is not the right way to study it -- the force fields are not that accurate. On Tue, Feb 25, 2014 at 6:48 AM, Dr. Vitaly Chaban vvcha...@gmail.com wrote: in that case the average is reasonable, but you should state in the manuscript that all the reported results correspond to 5.15 bar. Dr. Vitaly V. Chaban On Tue, Feb 25, 2014 at 10:42 AM, sujithkakkat . sujithk...@gmail.com wrote: Hi, In fact the average pressure 5.15 bar was obtained for the case where the reference pressure was 5 bar. I had mentioned it in one of the previous mails. But still the average value is not that good. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. Sujith. On Tue, Feb 25, 2014 at 2:29 PM, Dr. Vitaly Chaban vvcha...@gmail.com wrote: your average pressure is the pressure you should report in the publication. if you got 5 bars instead of 1 bar, you should write I simulated the system at 5 bars Dr. Vitaly V. Chaban On Tue, Feb 25, 2014 at 6:01 AM, sujithkakkat . sujithk...@gmail.com wrote: On Mon, Feb 24, 2014 at 5:37 AM, sujithkakkat . sujithk...@gmail.com wrote: Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. No, this average pressure is not fine. Why isn't it fine? I admit it is not very good. I guess due to huge fluctuations the average pressure is tend to vary a bit to either side of the reference value, unless I continue the equilibration for a very long time. Also the average was calculated with g_energy, which I guess considered the entire simulation, where the first few time steps, when the system is away from equilibrium, would have poor pressure values with respect to reference , which would reflect in the average pressure calculated. Also using g_analyze to calculate average pressure between different times during the simulation did not give the same value. So I am not sure if I should aim at an average pressure which is equal to the reference pressure. Thanks, Sujith. On Mon, Feb 24, 2014 at 4:26 PM, Justin Lemkul jalem...@vt.edu wrote: On 2/23/14, 11:37 PM, sujithkakkat . wrote: Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. I guess the sigma for the case of water was low and therefore the small cut-off of 0.8nm did not matter. However the case of cyclohexane alone remains to be tried. I guess Dr Vitaly was saying about using a switch/shift function. I will try the simulation with the new settings and see. The Gromos force fields were parametrized with a plain cutoff for van der Waals. They do not require switching or shifting. It's true that plain cutoffs are a bit harsh, but switching and shifting have their own issues, and one should always stick to the parametrization methods unless there is demonstrable evidence that the new method is better. I have never known Gromos96 force fields to need anything other than
Re: [gmx-users] Dubious results with NPT
On Tue, Feb 25, 2014 at 3:57 PM, Michael Shirts mrshi...@gmail.com wrote: No full solutions to this problem -- I'll write a few notes. * With large systems run for relatively small amounts of time, the ensemble could be statistically indistinguishble from the true distribution even if the averages don't line up correctly. See: http://pubs.acs.org/doi/abs/10.1021/ct300688p For ways to examine the full distribution of volumes. * Parrinello-Rahman isn't quite correct in GROMACS because virial and kinetic energy are offset by 1/2 a step, though getting something more like 1.2 atm when running at 1 atm is more common. I have been ever wondering why 1.2 and not 0.8 (from the own practice). * Note that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. * If you have a property that changes when the pressure changes from 1 to 5 atm, odds are that molecular simulation is not the right way to study it -- the force fields are not that accurate. I believe this is not a force field problem. We calculate normal boiling points by plotting saturated vapor pressure versus temperature. The results are quite trustworthy. It is NVT simulation, of course. As the problems with simulated pressure (and conjugated properties) are, to a large extent, due to a nanoscale size of the simulated systems, it is probably worthwhile to estimate system density in certain manner and then proceed with it in the NVT. Dr. Vitaly V. Chaban On Tue, Feb 25, 2014 at 6:48 AM, Dr. Vitaly Chaban vvcha...@gmail.com wrote: in that case the average is reasonable, but you should state in the manuscript that all the reported results correspond to 5.15 bar. Dr. Vitaly V. Chaban On Tue, Feb 25, 2014 at 10:42 AM, sujithkakkat . sujithk...@gmail.com wrote: Hi, In fact the average pressure 5.15 bar was obtained for the case where the reference pressure was 5 bar. I had mentioned it in one of the previous mails. But still the average value is not that good. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. Sujith. On Tue, Feb 25, 2014 at 2:29 PM, Dr. Vitaly Chaban vvcha...@gmail.com wrote: your average pressure is the pressure you should report in the publication. if you got 5 bars instead of 1 bar, you should write I simulated the system at 5 bars Dr. Vitaly V. Chaban On Tue, Feb 25, 2014 at 6:01 AM, sujithkakkat . sujithk...@gmail.com wrote: On Mon, Feb 24, 2014 at 5:37 AM, sujithkakkat . sujithk...@gmail.com wrote: Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. No, this average pressure is not fine. Why isn't it fine? I admit it is not very good. I guess due to huge fluctuations the average pressure is tend to vary a bit to either side of the reference value, unless I continue the equilibration for a very long time. Also the average was calculated with g_energy, which I guess considered the entire simulation, where the first few time steps, when the system is away from equilibrium, would have poor pressure values with respect to reference , which would reflect in the average pressure calculated. Also using g_analyze to calculate average pressure between different times during the simulation did not give the same value. So I am not sure if I should aim at an average pressure which is equal to the reference pressure. Thanks, Sujith. On Mon, Feb 24, 2014 at 4:26 PM, Justin Lemkul jalem...@vt.edu wrote: On 2/23/14, 11:37 PM, sujithkakkat . wrote: Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. I guess the sigma for the case of water was low and therefore the small cut-off of 0.8nm did not
Re: [gmx-users] Dubious results with NPT
On 2/23/14, 11:37 PM, sujithkakkat . wrote: Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. I guess the sigma for the case of water was low and therefore the small cut-off of 0.8nm did not matter. However the case of cyclohexane alone remains to be tried. I guess Dr Vitaly was saying about using a switch/shift function. I will try the simulation with the new settings and see. The Gromos force fields were parametrized with a plain cutoff for van der Waals. They do not require switching or shifting. It's true that plain cutoffs are a bit harsh, but switching and shifting have their own issues, and one should always stick to the parametrization methods unless there is demonstrable evidence that the new method is better. I have never known Gromos96 force fields to need anything other than plain cutoff for van der Waals to function as expected. -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 601 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Dubious results with NPT
On 2/23/14, 8:30 AM, sujith wrote: Dear GROMACS users, I am new to GROMACS, and recently started using the version 4.6.5. I have seen a lot of NPT related issues raised earlier in this forum, but in my case the error looks much more severe. I am following Justin Lemkul's tutorial (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/biphasic/01_genconf.html) on Biphasic system, containing cyclohexane and water.Everything went well till the NPT simulation to bring system to reference pressure of 1 bar. Here are the details of the system , .mdp file , what I did and where I find the problem. SYSTEM : 512 cyclohexane + 2720 water molecules. CURRENT STATUS: (1)Energy minimization : energy converged. (2)NVT (tcoupl=berendsen , tau_t=0.1ps / ref_t=298K) completed and target pressure of 298 K attained. (3)NPT (tcoupl=nose-hoover, tau_t=0.5ps / pcoupl=berendsen, tau_p=1ps / ref_p = 1bar , total time=15 ns) completed with the following result: Energy Average Err.Est. RMSD Tot-Drift --- Pressure1.08924 0.67 114.66 -1.74033 (bar) Since berendsen barostat doesn't generate the true ensemble, an equilibration with pcoupl=parrinello-rahman is performed in the next step. (4) NPT (pcoupl=parrinello-rahman, tau_p=5ps, total time=5ns) PROBLEM FACED: The average pressure too high as shown below which I feel is not going to improve. Energy Average Err.Est. RMSD Tot-Drift --- Pressure23.6651 0.53144.464 -1.00634 (bar) I am aware of the fact that pressure is subject to large fluctuations in small sized systems and that this may affect the average value of the pressure. But , here the average pressure looks too large to be ignored. The pressure-vs-time graph doesn't show any upward trend, and the pressure looks like fluctuating about a central value. Here is the .mdp file. Only changes made from the previous .mdp for berendsen pressure coupling , are with pcoupl and tau_p. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.002 nsteps = 250 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 2500 nstvout = 2500 nstfout = 2500 nstlog = 2500 nstenergy = 100 nstxtcout = 1000 xtc_precision = 1000 xtc_grps= System energygrps = System ; 7.3.9 Neighbor Searching nstlist = 1 ns_type= grid pbc = xyz rlist = 0.8 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.8 ; 7.3.11 VdW vdwtype = cut-off rvdw = 0.8 What force field are you using? If it's Gromos96 like my tutorial, the value of rvdw is much too short and can lead to artifacts. Using 0.8 for rlist/rcoulomb is fine, though 0.9 is more common. rvdw should be 1.4. DispCorr= EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol = 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t= 0.50.5 ref_t= 298298 ; 7.3.15 Pressure Coupling pcoupl= parrinello-rahman pcoupltype = isotropic tau_p = 5.0 compressibility= 4.5e-5 ref_p = 1.0 The value of tau_p seems a bit long to me; does changing it to 1.0 improve the results? ; 7.3.17 Velocity Generation gen_vel = no ; 7.3.18 Bonds constraints = all-bonds constraint_algorithm= LINCS continuation = yes lincs_order = 4 lincs_iter= 1 lincs_warnangle = 30 npt.mdp 63L, 4080C I guess there is something seriously wrong in the choice of methods/parameters in the .mdp file, which I cant figure out. Kindly go through and let me know your comments. I would be happy to give any further details. Any help would be appreciated. An even better test is to simplify the system. Run pure water or pure cyclohexane with the existing settings, then try the modifications above. That should help root out whether you're having a problem with the .mdp settings, force field, or maybe even the combination of both. -Justin -- == Justin A. Lemkul,
Re: [gmx-users] Dubious results with NPT
Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. I guess the sigma for the case of water was low and therefore the small cut-off of 0.8nm did not matter. However the case of cyclohexane alone remains to be tried. I guess Dr Vitaly was saying about using a switch/shift function. I will try the simulation with the new settings and see. Regards, Sujith. -- View this message in context: http://gromacs.5086.x6.nabble.com/Dubious-results-with-NPT-tp5014731p5014759.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Dubious results with NPT
Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. I guess the sigma for the case of water was low and therefore the small cut-off of 0.8nm did not matter. However the case of cyclohexane alone remains to be tried. I guess Dr Vitaly was saying about using a switch/shift function. I will try the simulation with the new settings and see. Regards, Sujith. On Sun, Feb 23, 2014 at 10:37 PM, Dr. Vitaly Chaban vvcha...@gmail.comwrote: There is such thing in simulations as energy conservation... If you use vdwtype= cut-off and this cut-off happens at 0.8nm, while sigma for the largest atom is ~0.34nm, the problems are inevitable. Your cut-off should not be smaller than 0.90nm, and you need to apply a more polite method to bring pairwise energy term down to zero at this cut-off. Dr. Vitaly V. Chaban On Sun, Feb 23, 2014 at 3:07 PM, Justin Lemkul jalem...@vt.edu wrote: On 2/23/14, 8:30 AM, sujith wrote: Dear GROMACS users, I am new to GROMACS, and recently started using the version 4.6.5. I have seen a lot of NPT related issues raised earlier in this forum, but in my case the error looks much more severe. I am following Justin Lemkul's tutorial ( http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/biphasic/01_genconf.html ) on Biphasic system, containing cyclohexane and water.Everything went well till the NPT simulation to bring system to reference pressure of 1 bar. Here are the details of the system , .mdp file , what I did and where I find the problem. SYSTEM : 512 cyclohexane + 2720 water molecules. CURRENT STATUS: (1)Energy minimization : energy converged. (2)NVT (tcoupl=berendsen , tau_t=0.1ps / ref_t=298K) completed and target pressure of 298 K attained. (3)NPT (tcoupl=nose-hoover, tau_t=0.5ps / pcoupl=berendsen, tau_p=1ps / ref_p = 1bar , total time=15 ns) completed with the following result: Energy Average Err.Est. RMSD Tot-Drift --- Pressure1.08924 0.67 114.66 -1.74033 (bar) Since berendsen barostat doesn't generate the true ensemble, an equilibration with pcoupl=parrinello-rahman is performed in the next step. (4) NPT (pcoupl=parrinello-rahman, tau_p=5ps, total time=5ns) PROBLEM FACED: The average pressure too high as shown below which I feel is not going to improve. Energy Average Err.Est. RMSD Tot-Drift --- Pressure23.6651 0.53144.464 -1.00634 (bar) I am aware of the fact that pressure is subject to large fluctuations in small sized systems and that this may affect the average value of the pressure. But , here the average pressure looks too large to be ignored. The pressure-vs-time graph doesn't show any upward trend, and the pressure looks like fluctuating about a central value. Here is the .mdp file. Only changes made from the previous .mdp for berendsen pressure coupling , are with pcoupl and tau_p. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.002 nsteps = 250 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 2500 nstvout = 2500 nstfout = 2500 nstlog = 2500 nstenergy = 100 nstxtcout = 1000 xtc_precision = 1000 xtc_grps= System energygrps = System ; 7.3.9 Neighbor Searching nstlist = 1 ns_type= grid pbc = xyz rlist = 0.8 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.8 ; 7.3.11 VdW vdwtype = cut-off rvdw = 0.8 What force field are you using? If it's Gromos96 like my tutorial, the value of rvdw is much too short and can lead to artifacts. Using 0.8 for
