Re: [gmx-users] Is anyone also using lammps?s
Dear Peng, Did you also try to run GMX in double precision at some point? Ran Peng Yi wrote: I turned off the torsion interaction. The difference between Lammps and Gromacs at integration time step 2fs was reduced. Details below: A melt of 240 n-octane (united-atom model), NVT, T=300K, V=55.46nm^3. Both Lammps and Gromacs use berendsen thermostat with tau_t=1ps. Integration time step 1fs: LammpsGromacsStd. Err. (for both) Ebond(kJ/mol):2092 2109 100 Eangle: 1778 175480 Elj+corr: -10501-10553 100 T(K): 300 299 5 P(atm): 3188 3016 700 integration time step 2fs: LammpsGromacsStd.Err. (for both) Ebond:2133 2232 100 Eangle: 1803 173780 Elj+corr: -10501-10623 100 T: 300 298 5 P:3133 2955 600 Lammps results remain almost unchanged when dt increases from 1fs to 2fs, and Ebond : Eangle = 7 : 6, which is the ratio between # of bonds and # of angles. Gromacs results change more significantly when dt goes from 1fs to 2fs. and the trend of Ebond and Eangle are opposite. It is more significant when torsion interaction is present. On Tue, 3 Nov 2009, David van der Spoel wrote: Peng Yi wrote: Hi, David, I used Berendsen thermostat and a bigger tau_t=1ps to redo the simulations. The general conclusion is the same. The std err is the same in both packages. And during each simulation, the integration time does not change. Details below: Integration time step 1fs: Lammps GromacsStd.Err. (for both) Ebond(kJ/mol): 21122170 100 Eangle: 17991770 100 Etors: 25522490 100 Elj+corr: -10711 -10777 100 P(atm): 32503216 500 Integration time step 2fs: Lammps GromacsStd.Err. (for both) Ebond: 21542654 120 Eangle: 18401645 120 Etors: 25732236 120 Elj+corr: -10711 -11019 100 P(atm): 3250 2590 600 How about the temperature in both systems? Was Lammps also run with Berendsen? It could also still be a topology error. Maybe you can turn off the torsion potential to test this. -Peng On Sun, 1 Nov 2009, David van der Spoel wrote: Peng Yi wrote: Thank for your reply! I have done some NVT runs per your suggestion, and the results are similar to NPT runs, i.e., Gromacs results is more affected by changing integration timestep than Lammps. Details below: A melt of 240 octane chains by united-atom model. T=300K, V=55.46 nm^3. Both Gromacs and Lammps use Nose-Hoover thermostat with tau_t=0.2 ps. As some people have note the tau_t is short for Nose Hoover. Are you sure this means the same in Lammps and in Gromacs? For one thing, there is no tau_t in the NH algorithm as far as I know, and Gromacs converts it to an appropriate weight or whatever that is called. What does Lammps do with this tau_t. To be a the safe side you could run both with Berendsen as well. Is the std err identical in both packages? And in the 2 fs run, are both simulation equlibrated with this time step as well? All other parameters in .top and .mdp files are the same as previously attached.. If I use integration time 1fs, Lammps and Gromacs produce consistent results: Lammps Gromacs std. err. Ebond(kJ/mol):21332160 100 Eangle: 17571780 80 Etors:25312510 80 Elj+corr: -10711 -10767 90 P(atm): 35003250 500 if I use integration time 2fs, Lammps results remain unchanged, but Gromacs results change significantly, particularly bonded energy: Lammps Gromacs std. err. Ebond(kJ/mol):21752710 100 Eangle: 17991640 70 Etors:25732230 80 Elj+corr: -10711 -11007 100 P(atm): 3200 2730 700 Would that be a result of using different integrator between Lammps and Gromacs? Lammps uses Velocity-Verlet, and Gromacs uses Leap-frog. Thanks, -Peng On Thu, 29 Oct 2009, David van der Spoel wrote: aherz wrote: Hey, are you running single or double precision gromacs? Afaik, depending on the circumstances the energy drift in gromacs can be rather bad for single precision. Please refer to the gromacs 4.0 paper for a discussion of the drift. If you want to compare energies you need the same density, which you do not have, you may need to run NVT for
Re: [gmx-users] Is anyone also using lammps?s
Hi, Ran, No, I haven't. I still have to find out how to install in double precision. Would double precision be slower than single? If so, how much? Or just double the memory used? Thanks, -Peng On Wed, 4 Nov 2009, Ran Friedman wrote: Dear Peng, Did you also try to run GMX in double precision at some point? Ran Peng Yi wrote: I turned off the torsion interaction. The difference between Lammps and Gromacs at integration time step 2fs was reduced. Details below: A melt of 240 n-octane (united-atom model), NVT, T=300K, V=55.46nm^3. Both Lammps and Gromacs use berendsen thermostat with tau_t=1ps. Integration time step 1fs: LammpsGromacsStd. Err. (for both) Ebond(kJ/mol):2092 2109 100 Eangle: 1778 175480 Elj+corr: -10501-10553 100 T(K): 300 299 5 P(atm): 3188 3016 700 integration time step 2fs: LammpsGromacsStd.Err. (for both) Ebond:2133 2232 100 Eangle: 1803 173780 Elj+corr: -10501-10623 100 T: 300 298 5 P:3133 2955 600 Lammps results remain almost unchanged when dt increases from 1fs to 2fs, and Ebond : Eangle = 7 : 6, which is the ratio between # of bonds and # of angles. Gromacs results change more significantly when dt goes from 1fs to 2fs. and the trend of Ebond and Eangle are opposite. It is more significant when torsion interaction is present. On Tue, 3 Nov 2009, David van der Spoel wrote: Peng Yi wrote: Hi, David, I used Berendsen thermostat and a bigger tau_t=1ps to redo the simulations. The general conclusion is the same. The std err is the same in both packages. And during each simulation, the integration time does not change. Details below: Integration time step 1fs: Lammps GromacsStd.Err. (for both) Ebond(kJ/mol): 21122170 100 Eangle: 17991770 100 Etors: 25522490 100 Elj+corr: -10711 -10777 100 P(atm): 32503216 500 Integration time step 2fs: Lammps GromacsStd.Err. (for both) Ebond: 21542654 120 Eangle: 18401645 120 Etors: 25732236 120 Elj+corr: -10711 -11019 100 P(atm): 3250 2590 600 How about the temperature in both systems? Was Lammps also run with Berendsen? It could also still be a topology error. Maybe you can turn off the torsion potential to test this. -Peng On Sun, 1 Nov 2009, David van der Spoel wrote: Peng Yi wrote: Thank for your reply! I have done some NVT runs per your suggestion, and the results are similar to NPT runs, i.e., Gromacs results is more affected by changing integration timestep than Lammps. Details below: A melt of 240 octane chains by united-atom model. T=300K, V=55.46 nm^3. Both Gromacs and Lammps use Nose-Hoover thermostat with tau_t=0.2 ps. As some people have note the tau_t is short for Nose Hoover. Are you sure this means the same in Lammps and in Gromacs? For one thing, there is no tau_t in the NH algorithm as far as I know, and Gromacs converts it to an appropriate weight or whatever that is called. What does Lammps do with this tau_t. To be a the safe side you could run both with Berendsen as well. Is the std err identical in both packages? And in the 2 fs run, are both simulation equlibrated with this time step as well? All other parameters in .top and .mdp files are the same as previously attached.. If I use integration time 1fs, Lammps and Gromacs produce consistent results: Lammps Gromacs std. err. Ebond(kJ/mol):21332160 100 Eangle: 17571780 80 Etors:25312510 80 Elj+corr: -10711 -10767 90 P(atm): 35003250 500 if I use integration time 2fs, Lammps results remain unchanged, but Gromacs results change significantly, particularly bonded energy: Lammps Gromacs std. err. Ebond(kJ/mol):21752710 100 Eangle: 17991640 70 Etors:25732230 80 Elj+corr: -10711 -11007 100 P(atm): 3200 2730 700 Would that be a result of using different integrator between Lammps and Gromacs? Lammps uses Velocity-Verlet, and Gromacs uses Leap-frog. Thanks, -Peng On Thu, 29 Oct 2009, David van der Spoel wrote: aherz wrote: Hey, are you running single or double precision gromacs? Afaik, depending on the circumstances the energy drift in gromacs can be rather bad for single precision. Please refer to the gromacs 4.0 paper
Re: [gmx-users] Is anyone also using lammps?s
Hi, Ran, Do you mean the simulation I used here for testing purpose? That is not long, maybe a few hours. But my research will require simulations for days. If I will not use PME, do I still need a fftw? -Peng On Wed, 4 Nov 2009, Ran Friedman wrote: Hi Peng, It would be slower - never made any benchmarks on how much slower. But you don't run a very long simulation, do you? Installing it isn't a problem. If you use PME you also need fftw in double precision though. Ran Peng Yi wrote: Hi, Ran, No, I haven't. I still have to find out how to install in double precision. Would double precision be slower than single? If so, how much? Or just double the memory used? Thanks, -Peng ___ 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] Is anyone also using lammps?s
Hi Peng, AFAIK GMX only uses fftw for PME. If you have no reason to prefer flexible bonds you can use LINCS or SHAKE to constrain the bond lengths and run with a timestep of 2fs at room temperature. You may want to try LAMMPS and GMX with SHAKE and a timestep of 2fs and see if the results are similar before recompiling GMX. If I understood you correctly you didn't constrain the bonds. Ran Peng Yi wrote: Hi, Ran, Do you mean the simulation I used here for testing purpose? That is not long, maybe a few hours. But my research will require simulations for days. If I will not use PME, do I still need a fftw? -Peng On Wed, 4 Nov 2009, Ran Friedman wrote: Hi Peng, It would be slower - never made any benchmarks on how much slower. But you don't run a very long simulation, do you? Installing it isn't a problem. If you use PME you also need fftw in double precision though. Ran Peng Yi wrote: Hi, Ran, No, I haven't. I still have to find out how to install in double precision. Would double precision be slower than single? If so, how much? Or just double the memory used? Thanks, -Peng ___ 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 ___ 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] Is anyone also using lammps?s
Hi, David, I used Berendsen thermostat and a bigger tau_t=1ps to redo the simulations. The general conclusion is the same. The std err is the same in both packages. And during each simulation, the integration time does not change. Details below: Integration time step 1fs: Lammps GromacsStd.Err. (for both) Ebond(kJ/mol): 21122170 100 Eangle: 17991770 100 Etors: 25522490 100 Elj+corr: -10711 -10777 100 P(atm): 32503216 500 Integration time step 2fs: Lammps GromacsStd.Err. (for both) Ebond: 21542654 120 Eangle: 18401645 120 Etors: 25732236 120 Elj+corr: -10711 -11019 100 P(atm): 32502590 600 -Peng On Sun, 1 Nov 2009, David van der Spoel wrote: Peng Yi wrote: Thank for your reply! I have done some NVT runs per your suggestion, and the results are similar to NPT runs, i.e., Gromacs results is more affected by changing integration timestep than Lammps. Details below: A melt of 240 octane chains by united-atom model. T=300K, V=55.46 nm^3. Both Gromacs and Lammps use Nose-Hoover thermostat with tau_t=0.2 ps. As some people have note the tau_t is short for Nose Hoover. Are you sure this means the same in Lammps and in Gromacs? For one thing, there is no tau_t in the NH algorithm as far as I know, and Gromacs converts it to an appropriate weight or whatever that is called. What does Lammps do with this tau_t. To be a the safe side you could run both with Berendsen as well. Is the std err identical in both packages? And in the 2 fs run, are both simulation equlibrated with this time step as well? All other parameters in .top and .mdp files are the same as previously attached.. If I use integration time 1fs, Lammps and Gromacs produce consistent results: Lammps Gromacs std. err. Ebond(kJ/mol):21332160 100 Eangle: 17571780 80 Etors:25312510 80 Elj+corr: -10711 -10767 90 P(atm): 35003250 500 if I use integration time 2fs, Lammps results remain unchanged, but Gromacs results change significantly, particularly bonded energy: Lammps Gromacs std. err. Ebond(kJ/mol):21752710 100 Eangle: 17991640 70 Etors:25732230 80 Elj+corr: -10711 -11007 100 P(atm): 32002730 700 Would that be a result of using different integrator between Lammps and Gromacs? Lammps uses Velocity-Verlet, and Gromacs uses Leap-frog. Thanks, -Peng On Thu, 29 Oct 2009, David van der Spoel wrote: aherz wrote: Hey, are you running single or double precision gromacs? Afaik, depending on the circumstances the energy drift in gromacs can be rather bad for single precision. Please refer to the gromacs 4.0 paper for a discussion of the drift. If you want to compare energies you need the same density, which you do not have, you may need to run NVT for that. Note that your integration time step is quite large, and the temperature coupling constant is very small. You could try a shifted LJ + dispersion correction, it is not clear to me how LAMMPS treats cutoffs, couldn't find it in the manual. Alex Peng Yi schrieb: On Wed, 28 Oct 2009, Mark Abraham wrote: Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark Hi, Mark, Thanks for reply! The difference is statistically significant. And I am wondering if it is caused by the integrator: Leap-frog for Gromacs and Velocity-verlet for Lammps. Detail description of the comparison please see below: It is an NPT simulation of a melt of 240 n-octane molecules using united-atom model, i.e., CHx group is considered as one atom. There are bond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs uses nose-hoover thermostat and Parranello-Rahman barostat. Time constants for thermostat and barostat are 0.02ps and 2.0ps, respectively. If I use integration time 1fs, Lammps and Gromacs gave consistent results: Lammps Gromacs Ebond(kJ/mol):2092 2146 Eangle: 1757 1760 Etors:2510 2500
Re: [gmx-users] Is anyone also using lammps?s
Peng Yi wrote: Hi, David, I used Berendsen thermostat and a bigger tau_t=1ps to redo the simulations. The general conclusion is the same. The std err is the same in both packages. And during each simulation, the integration time does not change. Details below: Integration time step 1fs: Lammps GromacsStd.Err. (for both) Ebond(kJ/mol): 21122170 100 Eangle: 17991770 100 Etors: 25522490 100 Elj+corr: -10711 -10777 100 P(atm): 32503216 500 Integration time step 2fs: Lammps GromacsStd.Err. (for both) Ebond: 21542654 120 Eangle: 18401645 120 Etors: 25732236 120 Elj+corr: -10711 -11019 100 P(atm): 32502590 600 How about the temperature in both systems? Was Lammps also run with Berendsen? It could also still be a topology error. Maybe you can turn off the torsion potential to test this. -Peng On Sun, 1 Nov 2009, David van der Spoel wrote: Peng Yi wrote: Thank for your reply! I have done some NVT runs per your suggestion, and the results are similar to NPT runs, i.e., Gromacs results is more affected by changing integration timestep than Lammps. Details below: A melt of 240 octane chains by united-atom model. T=300K, V=55.46 nm^3. Both Gromacs and Lammps use Nose-Hoover thermostat with tau_t=0.2 ps. As some people have note the tau_t is short for Nose Hoover. Are you sure this means the same in Lammps and in Gromacs? For one thing, there is no tau_t in the NH algorithm as far as I know, and Gromacs converts it to an appropriate weight or whatever that is called. What does Lammps do with this tau_t. To be a the safe side you could run both with Berendsen as well. Is the std err identical in both packages? And in the 2 fs run, are both simulation equlibrated with this time step as well? All other parameters in .top and .mdp files are the same as previously attached.. If I use integration time 1fs, Lammps and Gromacs produce consistent results: Lammps Gromacs std. err. Ebond(kJ/mol):21332160 100 Eangle: 17571780 80 Etors:25312510 80 Elj+corr: -10711 -10767 90 P(atm): 35003250 500 if I use integration time 2fs, Lammps results remain unchanged, but Gromacs results change significantly, particularly bonded energy: Lammps Gromacs std. err. Ebond(kJ/mol):21752710 100 Eangle: 17991640 70 Etors:25732230 80 Elj+corr: -10711 -11007 100 P(atm): 32002730 700 Would that be a result of using different integrator between Lammps and Gromacs? Lammps uses Velocity-Verlet, and Gromacs uses Leap-frog. Thanks, -Peng On Thu, 29 Oct 2009, David van der Spoel wrote: aherz wrote: Hey, are you running single or double precision gromacs? Afaik, depending on the circumstances the energy drift in gromacs can be rather bad for single precision. Please refer to the gromacs 4.0 paper for a discussion of the drift. If you want to compare energies you need the same density, which you do not have, you may need to run NVT for that. Note that your integration time step is quite large, and the temperature coupling constant is very small. You could try a shifted LJ + dispersion correction, it is not clear to me how LAMMPS treats cutoffs, couldn't find it in the manual. Alex Peng Yi schrieb: On Wed, 28 Oct 2009, Mark Abraham wrote: Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark Hi, Mark, Thanks for reply! The difference is statistically significant. And I am wondering if it is caused by the integrator: Leap-frog for Gromacs and Velocity-verlet for Lammps. Detail description of the comparison please see below: It is an NPT simulation of a melt of 240 n-octane molecules using united-atom model, i.e., CHx group is considered as one atom. There are bond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs uses nose-hoover thermostat and Parranello-Rahman barostat. Time constants for thermostat and barostat are 0.02ps and 2.0ps, respectively. If I use integration time 1fs, Lammps and
Re: [gmx-users] Is anyone also using lammps?s
I turned off the torsion interaction. The difference between Lammps and Gromacs at integration time step 2fs was reduced. Details below: A melt of 240 n-octane (united-atom model), NVT, T=300K, V=55.46nm^3. Both Lammps and Gromacs use berendsen thermostat with tau_t=1ps. Integration time step 1fs: LammpsGromacsStd. Err. (for both) Ebond(kJ/mol):2092 2109 100 Eangle: 1778 175480 Elj+corr: -10501-10553 100 T(K): 300 299 5 P(atm): 3188 3016 700 integration time step 2fs: LammpsGromacsStd.Err. (for both) Ebond:2133 2232 100 Eangle: 1803 173780 Elj+corr: -10501-10623 100 T: 300 298 5 P:3133 2955 600 Lammps results remain almost unchanged when dt increases from 1fs to 2fs, and Ebond : Eangle = 7 : 6, which is the ratio between # of bonds and # of angles. Gromacs results change more significantly when dt goes from 1fs to 2fs. and the trend of Ebond and Eangle are opposite. It is more significant when torsion interaction is present. On Tue, 3 Nov 2009, David van der Spoel wrote: Peng Yi wrote: Hi, David, I used Berendsen thermostat and a bigger tau_t=1ps to redo the simulations. The general conclusion is the same. The std err is the same in both packages. And during each simulation, the integration time does not change. Details below: Integration time step 1fs: Lammps GromacsStd.Err. (for both) Ebond(kJ/mol): 21122170 100 Eangle: 17991770 100 Etors: 25522490 100 Elj+corr: -10711 -10777 100 P(atm): 32503216 500 Integration time step 2fs: Lammps GromacsStd.Err. (for both) Ebond: 21542654 120 Eangle: 18401645 120 Etors: 25732236 120 Elj+corr: -10711 -11019 100 P(atm): 3250 2590 600 How about the temperature in both systems? Was Lammps also run with Berendsen? It could also still be a topology error. Maybe you can turn off the torsion potential to test this. -Peng On Sun, 1 Nov 2009, David van der Spoel wrote: Peng Yi wrote: Thank for your reply! I have done some NVT runs per your suggestion, and the results are similar to NPT runs, i.e., Gromacs results is more affected by changing integration timestep than Lammps. Details below: A melt of 240 octane chains by united-atom model. T=300K, V=55.46 nm^3. Both Gromacs and Lammps use Nose-Hoover thermostat with tau_t=0.2 ps. As some people have note the tau_t is short for Nose Hoover. Are you sure this means the same in Lammps and in Gromacs? For one thing, there is no tau_t in the NH algorithm as far as I know, and Gromacs converts it to an appropriate weight or whatever that is called. What does Lammps do with this tau_t. To be a the safe side you could run both with Berendsen as well. Is the std err identical in both packages? And in the 2 fs run, are both simulation equlibrated with this time step as well? All other parameters in .top and .mdp files are the same as previously attached.. If I use integration time 1fs, Lammps and Gromacs produce consistent results: Lammps Gromacs std. err. Ebond(kJ/mol):21332160 100 Eangle: 17571780 80 Etors:25312510 80 Elj+corr: -10711 -10767 90 P(atm): 35003250 500 if I use integration time 2fs, Lammps results remain unchanged, but Gromacs results change significantly, particularly bonded energy: Lammps Gromacs std. err. Ebond(kJ/mol):21752710 100 Eangle: 17991640 70 Etors:25732230 80 Elj+corr: -10711 -11007 100 P(atm): 3200 2730 700 Would that be a result of using different integrator between Lammps and Gromacs? Lammps uses Velocity-Verlet, and Gromacs uses Leap-frog. Thanks, -Peng On Thu, 29 Oct 2009, David van der Spoel wrote: aherz wrote: Hey, are you running single or double precision gromacs? Afaik, depending on the circumstances the energy drift in gromacs can be rather bad for single precision. Please refer to the gromacs 4.0 paper for a discussion of the drift. If you want to compare energies you need the same density, which you do not have, you may need to run NVT for that. Note that your integration time step is quite large, and the temperature coupling constant is very small. You could try a shifted LJ + dispersion
Re: [gmx-users] Is anyone also using lammps?s
Thank for your reply! I have done some NVT runs per your suggestion, and the results are similar to NPT runs, i.e., Gromacs results is more affected by changing integration timestep than Lammps. Details below: A melt of 240 octane chains by united-atom model. T=300K, V=55.46 nm^3. Both Gromacs and Lammps use Nose-Hoover thermostat with tau_t=0.2 ps. All other parameters in .top and .mdp files are the same as previously attached.. If I use integration time 1fs, Lammps and Gromacs produce consistent results: Lammps Gromacs std. err. Ebond(kJ/mol):21332160 100 Eangle: 17571780 80 Etors:25312510 80 Elj+corr: -10711 -10767 90 P(atm): 35003250 500 if I use integration time 2fs, Lammps results remain unchanged, but Gromacs results change significantly, particularly bonded energy: Lammps Gromacs std. err. Ebond(kJ/mol):21752710 100 Eangle: 17991640 70 Etors:25732230 80 Elj+corr: -10711 -11007 100 P(atm): 32002730 700 Would that be a result of using different integrator between Lammps and Gromacs? Lammps uses Velocity-Verlet, and Gromacs uses Leap-frog. Thanks, -Peng On Thu, 29 Oct 2009, David van der Spoel wrote: aherz wrote: Hey, are you running single or double precision gromacs? Afaik, depending on the circumstances the energy drift in gromacs can be rather bad for single precision. Please refer to the gromacs 4.0 paper for a discussion of the drift. If you want to compare energies you need the same density, which you do not have, you may need to run NVT for that. Note that your integration time step is quite large, and the temperature coupling constant is very small. You could try a shifted LJ + dispersion correction, it is not clear to me how LAMMPS treats cutoffs, couldn't find it in the manual. Alex Peng Yi schrieb: On Wed, 28 Oct 2009, Mark Abraham wrote: Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark Hi, Mark, Thanks for reply! The difference is statistically significant. And I am wondering if it is caused by the integrator: Leap-frog for Gromacs and Velocity-verlet for Lammps. Detail description of the comparison please see below: It is an NPT simulation of a melt of 240 n-octane molecules using united-atom model, i.e., CHx group is considered as one atom. There are bond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs uses nose-hoover thermostat and Parranello-Rahman barostat. Time constants for thermostat and barostat are 0.02ps and 2.0ps, respectively. If I use integration time 1fs, Lammps and Gromacs gave consistent results: Lammps Gromacs Ebond(kJ/mol):2092 2146 Eangle: 1757 1760 Etors:2510 2500 Elj+corr:-9238-9350 Volume(nm^3): 66.7 66.5 where energy fluctuation is 100 kJ/mol and volume fluctuation is 1 nm^3, Elj+corr is the total LJ energy including tail correction. However, if I use integration time 2fs, Lammps results do not change much, but Gromacs results changed a lot: Lammps Gromacs Ebond(kJ/mol):2133 2700 Eangle: 1799 1640 Etors:2552 2200 Elj+corr:-9292-9886 Volume: 66.7 64.0 The results given by Lammps is more reasonable because the Ebond should be equal to the total # of bonds times 1/2k_BT and Eangle should be equal to the total # of angles times 1/2k_BT. At T=300K, 1/2k_BT=1.25 kJ/mol. 240 n-octanes have total 1680 bonds and 1440 angles. The bond and angle interactions are both harmonic functions. Bond interaction constant kl=292880 kJ/mol/nm^2, corresponding to a bond ossilation period 16 fs. Is there something related to the integrator? Here I attached my grompp.mdp and topol.top files. ## grompp.mdp ## ; VARIOUS PREPROCESSING OPTIONS title= Yo cpp = /usr/bin/cpp include = define = ; RUN CONTROL PARAMETERS integrator = md tinit= 0 dt = 0.001 nsteps = 200 init_step= 0
Re: [gmx-users] Is anyone also using lammps?s
On Wed, 28 Oct 2009, Mark Abraham wrote: Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark Hi, Mark, Thanks for reply! The difference is statistically significant. And I am wondering if it is caused by the integrator: Leap-frog for Gromacs and Velocity-verlet for Lammps. Detail description of the comparison please see below: It is an NPT simulation of a melt of 240 n-octane molecules using united-atom model, i.e., CHx group is considered as one atom. There are bond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs uses nose-hoover thermostat and Parranello-Rahman barostat. Time constants for thermostat and barostat are 0.02ps and 2.0ps, respectively. If I use integration time 1fs, Lammps and Gromacs gave consistent results: Lammps Gromacs Ebond(kJ/mol):2092 2146 Eangle: 1757 1760 Etors:2510 2500 Elj+corr:-9238-9350 Volume(nm^3): 66.7 66.5 where energy fluctuation is 100 kJ/mol and volume fluctuation is 1 nm^3, Elj+corr is the total LJ energy including tail correction. However, if I use integration time 2fs, Lammps results do not change much, but Gromacs results changed a lot: Lammps Gromacs Ebond(kJ/mol):2133 2700 Eangle: 1799 1640 Etors:2552 2200 Elj+corr:-9292-9886 Volume: 66.7 64.0 The results given by Lammps is more reasonable because the Ebond should be equal to the total # of bonds times 1/2k_BT and Eangle should be equal to the total # of angles times 1/2k_BT. At T=300K, 1/2k_BT=1.25 kJ/mol. 240 n-octanes have total 1680 bonds and 1440 angles. The bond and angle interactions are both harmonic functions. Bond interaction constant kl=292880 kJ/mol/nm^2, corresponding to a bond ossilation period 16 fs. Is there something related to the integrator? Here I attached my grompp.mdp and topol.top files. ## grompp.mdp ## ; VARIOUS PREPROCESSING OPTIONS title= Yo cpp = /usr/bin/cpp include = define = ; RUN CONTROL PARAMETERS integrator = md tinit= 0 dt = 0.001 nsteps = 200 init_step= 0 comm-mode= Linear nstcomm = 1 comm-grps= ; OUTPUT CONTROL OPTIONS nstxout = 5000 nstvout = 5000 nstfout = 5000 nstcheckpoint= 1 nstlog = 1000 nstenergy= 1000 nstxtcout= 5000 xtc-precision= 1000 xtc-grps = energygrps = ; NEIGHBORSEARCHING PARAMETERS nstlist = 10 ns_type = grid pbc = xyz rlist= 1.0025 domain-decomposition = no ; OPTIONS FOR ELECTROSTATICS AND VDW coulombtype = Cut-off rcoulomb-switch = 0 rcoulomb = 1.0025 epsilon-r= 1 vdw-type = Cut-off rvdw-switch = 0 ; default rvdw = 1.0025 ; default 1 nm DispCorr = EnerPres ;table-extension = 1.5 fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order= 4 ewald_rtol = 1e-05 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS Tcoupl = nose-hoover tc-grps = System tau_t= 0.02 ref_t= 300.0 Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p= 2.0 compressibility = 4.5e-5 ref_p= 1.0 andersen_seed= 815131 ; GENERATE VELOCITIES FOR STARTUP RUN gen_vel = yes gen_temp = 300 gen_seed = 2009 ; OPTIONS FOR BONDS constraints = none constraint-algorithm = Lincs unconstrained-start = no Shake-SOR= no shake-tol= 1e-04 lincs-order = 4 lincs-iter = 1 lincs-warnangle = 30 morse= no ; ENERGY GROUP EXCLUSIONS ; Pairs of energy groups for which all
Re: [gmx-users] Is anyone also using lammps?s
Hey, are you running single or double precision gromacs? Afaik, depending on the circumstances the energy drift in gromacs can be rather bad for single precision. Alex Peng Yi schrieb: On Wed, 28 Oct 2009, Mark Abraham wrote: Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark Hi, Mark, Thanks for reply! The difference is statistically significant. And I am wondering if it is caused by the integrator: Leap-frog for Gromacs and Velocity-verlet for Lammps. Detail description of the comparison please see below: It is an NPT simulation of a melt of 240 n-octane molecules using united-atom model, i.e., CHx group is considered as one atom. There are bond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs uses nose-hoover thermostat and Parranello-Rahman barostat. Time constants for thermostat and barostat are 0.02ps and 2.0ps, respectively. If I use integration time 1fs, Lammps and Gromacs gave consistent results: Lammps Gromacs Ebond(kJ/mol):2092 2146 Eangle: 1757 1760 Etors:2510 2500 Elj+corr:-9238-9350 Volume(nm^3): 66.7 66.5 where energy fluctuation is 100 kJ/mol and volume fluctuation is 1 nm^3, Elj+corr is the total LJ energy including tail correction. However, if I use integration time 2fs, Lammps results do not change much, but Gromacs results changed a lot: Lammps Gromacs Ebond(kJ/mol):2133 2700 Eangle: 1799 1640 Etors:2552 2200 Elj+corr:-9292-9886 Volume: 66.7 64.0 The results given by Lammps is more reasonable because the Ebond should be equal to the total # of bonds times 1/2k_BT and Eangle should be equal to the total # of angles times 1/2k_BT. At T=300K, 1/2k_BT=1.25 kJ/mol. 240 n-octanes have total 1680 bonds and 1440 angles. The bond and angle interactions are both harmonic functions. Bond interaction constant kl=292880 kJ/mol/nm^2, corresponding to a bond ossilation period 16 fs. Is there something related to the integrator? Here I attached my grompp.mdp and topol.top files. ## grompp.mdp ## ; VARIOUS PREPROCESSING OPTIONS title= Yo cpp = /usr/bin/cpp include = define = ; RUN CONTROL PARAMETERS integrator = md tinit= 0 dt = 0.001 nsteps = 200 init_step= 0 comm-mode= Linear nstcomm = 1 comm-grps= ; OUTPUT CONTROL OPTIONS nstxout = 5000 nstvout = 5000 nstfout = 5000 nstcheckpoint= 1 nstlog = 1000 nstenergy= 1000 nstxtcout= 5000 xtc-precision= 1000 xtc-grps = energygrps = ; NEIGHBORSEARCHING PARAMETERS nstlist = 10 ns_type = grid pbc = xyz rlist= 1.0025 domain-decomposition = no ; OPTIONS FOR ELECTROSTATICS AND VDW coulombtype = Cut-off rcoulomb-switch = 0 rcoulomb = 1.0025 epsilon-r= 1 vdw-type = Cut-off rvdw-switch = 0; default rvdw = 1.0025; default 1 nm DispCorr = EnerPres ;table-extension = 1.5 fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order= 4 ewald_rtol = 1e-05 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS Tcoupl = nose-hoover tc-grps = System tau_t= 0.02 ref_t= 300.0 Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p= 2.0 compressibility = 4.5e-5 ref_p= 1.0 andersen_seed= 815131 ; GENERATE VELOCITIES FOR STARTUP RUN gen_vel = yes gen_temp = 300 gen_seed = 2009 ; OPTIONS FOR BONDS constraints = none constraint-algorithm
Re: [gmx-users] Is anyone also using lammps?s
aherz wrote: Hey, are you running single or double precision gromacs? Afaik, depending on the circumstances the energy drift in gromacs can be rather bad for single precision. Please refer to the gromacs 4.0 paper for a discussion of the drift. If you want to compare energies you need the same density, which you do not have, you may need to run NVT for that. Note that your integration time step is quite large, and the temperature coupling constant is very small. You could try a shifted LJ + dispersion correction, it is not clear to me how LAMMPS treats cutoffs, couldn't find it in the manual. Alex Peng Yi schrieb: On Wed, 28 Oct 2009, Mark Abraham wrote: Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark Hi, Mark, Thanks for reply! The difference is statistically significant. And I am wondering if it is caused by the integrator: Leap-frog for Gromacs and Velocity-verlet for Lammps. Detail description of the comparison please see below: It is an NPT simulation of a melt of 240 n-octane molecules using united-atom model, i.e., CHx group is considered as one atom. There are bond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs uses nose-hoover thermostat and Parranello-Rahman barostat. Time constants for thermostat and barostat are 0.02ps and 2.0ps, respectively. If I use integration time 1fs, Lammps and Gromacs gave consistent results: Lammps Gromacs Ebond(kJ/mol):2092 2146 Eangle: 1757 1760 Etors:2510 2500 Elj+corr:-9238-9350 Volume(nm^3): 66.7 66.5 where energy fluctuation is 100 kJ/mol and volume fluctuation is 1 nm^3, Elj+corr is the total LJ energy including tail correction. However, if I use integration time 2fs, Lammps results do not change much, but Gromacs results changed a lot: Lammps Gromacs Ebond(kJ/mol):2133 2700 Eangle: 1799 1640 Etors:2552 2200 Elj+corr:-9292-9886 Volume: 66.7 64.0 The results given by Lammps is more reasonable because the Ebond should be equal to the total # of bonds times 1/2k_BT and Eangle should be equal to the total # of angles times 1/2k_BT. At T=300K, 1/2k_BT=1.25 kJ/mol. 240 n-octanes have total 1680 bonds and 1440 angles. The bond and angle interactions are both harmonic functions. Bond interaction constant kl=292880 kJ/mol/nm^2, corresponding to a bond ossilation period 16 fs. Is there something related to the integrator? Here I attached my grompp.mdp and topol.top files. ## grompp.mdp ## ; VARIOUS PREPROCESSING OPTIONS title= Yo cpp = /usr/bin/cpp include = define = ; RUN CONTROL PARAMETERS integrator = md tinit= 0 dt = 0.001 nsteps = 200 init_step= 0 comm-mode= Linear nstcomm = 1 comm-grps= ; OUTPUT CONTROL OPTIONS nstxout = 5000 nstvout = 5000 nstfout = 5000 nstcheckpoint= 1 nstlog = 1000 nstenergy= 1000 nstxtcout= 5000 xtc-precision= 1000 xtc-grps = energygrps = ; NEIGHBORSEARCHING PARAMETERS nstlist = 10 ns_type = grid pbc = xyz rlist= 1.0025 domain-decomposition = no ; OPTIONS FOR ELECTROSTATICS AND VDW coulombtype = Cut-off rcoulomb-switch = 0 rcoulomb = 1.0025 epsilon-r= 1 vdw-type = Cut-off rvdw-switch = 0; default rvdw = 1.0025; default 1 nm DispCorr = EnerPres ;table-extension = 1.5 fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order= 4 ewald_rtol = 1e-05 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS Tcoupl = nose-hoover tc-grps = System tau_t= 0.02 ref_t= 300.0 Pcoupl = Parrinello-Rahman Pcoupltype = isotropic
Re: [gmx-users] Is anyone also using lammps?s
Hi Peng, Note that you're not using any bond constraints in Gromacs and a timestep of 2fs may be too long. Also, tau_t=0.02 seems too short for me. With 1fs timescale the agreement seem good enough, but you didn't include estimated errors so it's hard to tell. Also, I assume you run GMX in single and LAMMPS in double precision. Did you check for convergence? Ran Peng Yi wrote: On Wed, 28 Oct 2009, Mark Abraham wrote: Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark Hi, Mark, Thanks for reply! The difference is statistically significant. And I am wondering if it is caused by the integrator: Leap-frog for Gromacs and Velocity-verlet for Lammps. Detail description of the comparison please see below: It is an NPT simulation of a melt of 240 n-octane molecules using united-atom model, i.e., CHx group is considered as one atom. There are bond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs uses nose-hoover thermostat and Parranello-Rahman barostat. Time constants for thermostat and barostat are 0.02ps and 2.0ps, respectively. If I use integration time 1fs, Lammps and Gromacs gave consistent results: Lammps Gromacs Ebond(kJ/mol):2092 2146 Eangle: 1757 1760 Etors:2510 2500 Elj+corr:-9238-9350 Volume(nm^3): 66.7 66.5 where energy fluctuation is 100 kJ/mol and volume fluctuation is 1 nm^3, Elj+corr is the total LJ energy including tail correction. However, if I use integration time 2fs, Lammps results do not change much, but Gromacs results changed a lot: Lammps Gromacs Ebond(kJ/mol):2133 2700 Eangle: 1799 1640 Etors:2552 2200 Elj+corr:-9292-9886 Volume: 66.7 64.0 The results given by Lammps is more reasonable because the Ebond should be equal to the total # of bonds times 1/2k_BT and Eangle should be equal to the total # of angles times 1/2k_BT. At T=300K, 1/2k_BT=1.25 kJ/mol. 240 n-octanes have total 1680 bonds and 1440 angles. The bond and angle interactions are both harmonic functions. Bond interaction constant kl=292880 kJ/mol/nm^2, corresponding to a bond ossilation period 16 fs. Is there something related to the integrator? Here I attached my grompp.mdp and topol.top files. ## grompp.mdp ## ; VARIOUS PREPROCESSING OPTIONS title= Yo cpp = /usr/bin/cpp include = define = ; RUN CONTROL PARAMETERS integrator = md tinit= 0 dt = 0.001 nsteps = 200 init_step= 0 comm-mode= Linear nstcomm = 1 comm-grps= ; OUTPUT CONTROL OPTIONS nstxout = 5000 nstvout = 5000 nstfout = 5000 nstcheckpoint= 1 nstlog = 1000 nstenergy= 1000 nstxtcout= 5000 xtc-precision= 1000 xtc-grps = energygrps = ; NEIGHBORSEARCHING PARAMETERS nstlist = 10 ns_type = grid pbc = xyz rlist= 1.0025 domain-decomposition = no ; OPTIONS FOR ELECTROSTATICS AND VDW coulombtype = Cut-off rcoulomb-switch = 0 rcoulomb = 1.0025 epsilon-r= 1 vdw-type = Cut-off rvdw-switch = 0; default rvdw = 1.0025; default 1 nm DispCorr = EnerPres ;table-extension = 1.5 fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order= 4 ewald_rtol = 1e-05 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS Tcoupl = nose-hoover tc-grps = System tau_t= 0.02 ref_t= 300.0 Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p= 2.0 compressibility = 4.5e-5 ref_p= 1.0 andersen_seed= 815131 ; GENERATE
Re: [gmx-users] Is anyone also using lammps?s
Hi, Ran, I didn't use bond restraints. I checked that the bond length had a Gaussian-like distributes, and the length range looked normal. I estimated the fastest timescale in the system, which is the bond ossilation period, around 16fs. Would that require as integration timestep much smaller than 1fs? With the parameters I have, could you recommend a set of tau_t and tau_p? I did mention the fluctuation, 100 kJ/mol for energy and 1 nm^3 for volume. And I ran GMX in single. Not sure about Lammps, should be double. All measured physical quantities converged well. Would you expect differece if I compile GMX in double? Would that be much slower? -Peng On Thu, 29 Oct 2009, Ran Friedman wrote: Hi Peng, Note that you're not using any bond constraints in Gromacs and a timestep of 2fs may be too long. Also, tau_t=0.02 seems too short for me. With 1fs timescale the agreement seem good enough, but you didn't include estimated errors so it's hard to tell. Also, I assume you run GMX in single and LAMMPS in double precision. Did you check for convergence? Ran Peng Yi wrote: On Wed, 28 Oct 2009, Mark Abraham wrote: Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark Hi, Mark, Thanks for reply! The difference is statistically significant. And I am wondering if it is caused by the integrator: Leap-frog for Gromacs and Velocity-verlet for Lammps. Detail description of the comparison please see below: It is an NPT simulation of a melt of 240 n-octane molecules using united-atom model, i.e., CHx group is considered as one atom. There are bond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs uses nose-hoover thermostat and Parranello-Rahman barostat. Time constants for thermostat and barostat are 0.02ps and 2.0ps, respectively. If I use integration time 1fs, Lammps and Gromacs gave consistent results: Lammps Gromacs Ebond(kJ/mol):2092 2146 Eangle: 1757 1760 Etors:2510 2500 Elj+corr:-9238-9350 Volume(nm^3): 66.7 66.5 where energy fluctuation is 100 kJ/mol and volume fluctuation is 1 nm^3, Elj+corr is the total LJ energy including tail correction. However, if I use integration time 2fs, Lammps results do not change much, but Gromacs results changed a lot: Lammps Gromacs Ebond(kJ/mol):2133 2700 Eangle: 1799 1640 Etors:2552 2200 Elj+corr:-9292-9886 Volume: 66.7 64.0 The results given by Lammps is more reasonable because the Ebond should be equal to the total # of bonds times 1/2k_BT and Eangle should be equal to the total # of angles times 1/2k_BT. At T=300K, 1/2k_BT=1.25 kJ/mol. 240 n-octanes have total 1680 bonds and 1440 angles. The bond and angle interactions are both harmonic functions. Bond interaction constant kl=292880 kJ/mol/nm^2, corresponding to a bond ossilation period 16 fs. Is there something related to the integrator? Here I attached my grompp.mdp and topol.top files. ## grompp.mdp ## ; VARIOUS PREPROCESSING OPTIONS title= Yo cpp = /usr/bin/cpp include = define = ; RUN CONTROL PARAMETERS integrator = md tinit= 0 dt = 0.001 nsteps = 200 init_step= 0 comm-mode= Linear nstcomm = 1 comm-grps= ; OUTPUT CONTROL OPTIONS nstxout = 5000 nstvout = 5000 nstfout = 5000 nstcheckpoint= 1 nstlog = 1000 nstenergy= 1000 nstxtcout= 5000 xtc-precision= 1000 xtc-grps = energygrps = ; NEIGHBORSEARCHING PARAMETERS nstlist = 10 ns_type = grid pbc = xyz rlist= 1.0025 domain-decomposition = no ; OPTIONS FOR ELECTROSTATICS AND VDW coulombtype = Cut-off rcoulomb-switch = 0 rcoulomb = 1.0025 epsilon-r= 1 vdw-type = Cut-off rvdw-switch = 0; default rvdw = 1.0025; default 1 nm DispCorr = EnerPres ;table-extension = 1.5 fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0
Re: [gmx-users] Is anyone also using lammps?s
On Thu, 29 Oct 2009, David van der Spoel wrote: aherz wrote: Hey, are you running single or double precision gromacs? Afaik, depending on the circumstances the energy drift in gromacs can be rather bad for single precision. Please refer to the gromacs 4.0 paper for a discussion of the drift. If you want to compare energies you need the same density, which you do not have, you may need to run NVT for that. Note that your integration time step is quite large, and the temperature coupling constant is very small. You could try a shifted LJ + dispersion correction, it is not clear to me how LAMMPS treats cutoffs, couldn't find it in the manual. I tried NPT simulation with only LJ interaction. With tau_t=0.2ps for both Gromacs and Lammps, all the other parameters same as I mentioned before. At integration time 2fs, Gromacs and Lammps produce same results and agree with published data. (N=4000, T=2.0, P=0.5 - rou=0.3). That's why I focused my attention to the difference on bonded energy. I think I used a simple cutoff plus tail correction for energy and pressure for both Gromacs and Lammps. Alex Peng Yi schrieb: On Wed, 28 Oct 2009, Mark Abraham wrote: Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark Hi, Mark, Thanks for reply! The difference is statistically significant. And I am wondering if it is caused by the integrator: Leap-frog for Gromacs and Velocity-verlet for Lammps. Detail description of the comparison please see below: It is an NPT simulation of a melt of 240 n-octane molecules using united-atom model, i.e., CHx group is considered as one atom. There are bond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs uses nose-hoover thermostat and Parranello-Rahman barostat. Time constants for thermostat and barostat are 0.02ps and 2.0ps, respectively. If I use integration time 1fs, Lammps and Gromacs gave consistent results: Lammps Gromacs Ebond(kJ/mol):2092 2146 Eangle: 1757 1760 Etors:2510 2500 Elj+corr:-9238-9350 Volume(nm^3): 66.7 66.5 where energy fluctuation is 100 kJ/mol and volume fluctuation is 1 nm^3, Elj+corr is the total LJ energy including tail correction. However, if I use integration time 2fs, Lammps results do not change much, but Gromacs results changed a lot: Lammps Gromacs Ebond(kJ/mol):2133 2700 Eangle: 1799 1640 Etors:2552 2200 Elj+corr:-9292-9886 Volume: 66.7 64.0 The results given by Lammps is more reasonable because the Ebond should be equal to the total # of bonds times 1/2k_BT and Eangle should be equal to the total # of angles times 1/2k_BT. At T=300K, 1/2k_BT=1.25 kJ/mol. 240 n-octanes have total 1680 bonds and 1440 angles. The bond and angle interactions are both harmonic functions. Bond interaction constant kl=292880 kJ/mol/nm^2, corresponding to a bond ossilation period 16 fs. Is there something related to the integrator? Here I attached my grompp.mdp and topol.top files. ## grompp.mdp ## ; VARIOUS PREPROCESSING OPTIONS title= Yo cpp = /usr/bin/cpp include = define = ; RUN CONTROL PARAMETERS integrator = md tinit= 0 dt = 0.001 nsteps = 200 init_step= 0 comm-mode= Linear nstcomm = 1 comm-grps= ; OUTPUT CONTROL OPTIONS nstxout = 5000 nstvout = 5000 nstfout = 5000 nstcheckpoint= 1 nstlog = 1000 nstenergy= 1000 nstxtcout= 5000 xtc-precision= 1000 xtc-grps = energygrps = ; NEIGHBORSEARCHING PARAMETERS nstlist = 10 ns_type = grid pbc = xyz rlist= 1.0025 domain-decomposition = no ; OPTIONS FOR ELECTROSTATICS AND VDW coulombtype = Cut-off rcoulomb-switch = 0 rcoulomb = 1.0025 epsilon-r= 1 vdw-type = Cut-off rvdw-switch = 0; default rvdw= 1.0025; default 1 nm DispCorr = EnerPres ;table-extension = 1.5 fourierspacing = 0.12
Re: [gmx-users] Is anyone also using lammps?s
Hi Peng,
The time scale should be much shorter than the fastest vibration. A rule
of thumb from the reference below is a factor of ten, but it would
depend on the precision. Running with double precision is shorter but I
didn't make benchmarks (perhaps other users have).
Appropriate values of tau_t and tau_p have also been discussed in this
list (search for references by Berk). I tend to use something like
tau_t=0.2 and tau_p=1.0.
I advise you to follow David's suggestion as well and run with NVT.
Ran.
Reference:
@book{Becker2001,
Author = {Becker, O. M. and MacKerell, A. D. Jr. and Roux, B. and
Watanabe, M.},
Title = {Computational biochemistry and biophysics},
Publisher = {Dekker, M.},
Address = {New York},
Year = {2001} }
Peng Yi wrote:
Hi, Ran,
I didn't use bond restraints. I checked that the bond length had a
Gaussian-like distributes, and the length range looked normal.
I estimated the fastest timescale in the system, which is the bond
ossilation period, around 16fs. Would that require as integration
timestep much smaller than 1fs?
With the parameters I have, could you recommend a set of tau_t and tau_p?
I did mention the fluctuation, 100 kJ/mol for energy and 1 nm^3 for
volume. And I ran GMX in single. Not sure about Lammps, should be
double. All measured physical quantities converged well. Would you
expect differece if I compile GMX in double? Would that be much slower?
-Peng
On Thu, 29 Oct 2009, Ran Friedman wrote:
Hi Peng,
Note that you're not using any bond constraints in Gromacs and a
timestep of 2fs may be too long.
Also, tau_t=0.02 seems too short for me.
With 1fs timescale the agreement seem good enough, but you didn't
include estimated errors so it's hard to tell. Also, I assume you run
GMX in single and LAMMPS in double precision. Did you check for
convergence?
Ran
Peng Yi wrote:
On Wed, 28 Oct 2009, Mark Abraham wrote:
Peng Yi wrote:
I am trying to simulate alkane melt and found out that gromacs and
lammps gave different results, particularly the bonded interaction
energy.
I wonder if anyone has such experience. Thanks,
Even two installations of the same version of GROMACS can give
different results. The question is whether when using comparable
model physics you observe the same ensemble averages.
Mark
Hi, Mark,
Thanks for reply! The difference is statistically significant. And
I am
wondering if it is caused by the integrator: Leap-frog for Gromacs and
Velocity-verlet for Lammps. Detail description of the comparison
please
see below:
It is an NPT simulation of a melt of 240 n-octane molecules using
united-atom model, i.e., CHx group is considered as one atom. There
are
bond, angle, torsion and LJ interactions. T=300K and P=1atm.
Lammps uses nose-hoover thermostat and barostat, and Gromacs uses
nose-hoover thermostat and Parranello-Rahman barostat. Time constants
for
thermostat and barostat are 0.02ps and 2.0ps, respectively.
If I use integration time 1fs, Lammps and Gromacs gave consistent
results:
Lammps Gromacs
Ebond(kJ/mol):2092 2146
Eangle: 1757 1760
Etors:2510 2500
Elj+corr:-9238-9350
Volume(nm^3): 66.7 66.5
where energy fluctuation is 100 kJ/mol and volume fluctuation is 1
nm^3,
Elj+corr is the total LJ energy including tail correction.
However, if I use integration time 2fs, Lammps results do not change
much, but Gromacs results changed a lot:
Lammps Gromacs
Ebond(kJ/mol):2133 2700 Eangle:
1799 1640
Etors:2552 2200
Elj+corr:-9292-9886 Volume:
66.7 64.0
The results given by Lammps is more reasonable because the Ebond should
be equal to the total # of bonds times 1/2k_BT and Eangle should be
equal
to the total # of angles times 1/2k_BT. At T=300K, 1/2k_BT=1.25
kJ/mol.
240 n-octanes have total 1680 bonds and 1440 angles.
The bond and angle interactions are both harmonic functions. Bond
interaction constant kl=292880 kJ/mol/nm^2, corresponding to a bond
ossilation period 16 fs.
Is there something related to the integrator?
Here I attached my grompp.mdp and topol.top files.
##
grompp.mdp
##
; VARIOUS PREPROCESSING OPTIONS
title= Yo
cpp = /usr/bin/cpp
include = define =
; RUN CONTROL PARAMETERS
integrator = md
tinit= 0
dt = 0.001
nsteps = 200
init_step= 0
comm-mode= Linear
nstcomm = 1
comm-grps=
; OUTPUT CONTROL OPTIONS
nstxout = 5000
nstvout = 5000
nstfout = 5000
[gmx-users] Is anyone also using lammps?s
I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, ___ 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] Is anyone also using lammps?s
Peng Yi wrote: I am trying to simulate alkane melt and found out that gromacs and lammps gave different results, particularly the bonded interaction energy. I wonder if anyone has such experience. Thanks, Even two installations of the same version of GROMACS can give different results. The question is whether when using comparable model physics you observe the same ensemble averages. Mark ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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
