GROMACS Langevin integrator seems less stable than that of LAMMPS or HOOMD-BLUE

User discussions
1

GROMACS version: 2022.2
GROMACS modification: No
Dear Community,

I’m trying to simulate a polymer melt in reduce units using GROMACS. The system consists of 400 100-bead polymers with a bead density 0.85 \sigma^{-3}. The system has harmonic bond, harmonic angle and pure repulsive LJ (WCA) potential. The simulation details were from this paper https://pubs.acs.org/doi/10.1021/acs.macromol.5c03613?fig=tgr1&ref=pdf

in which the drag coefficient is 0.5 \tau^{-1} Thus, I set mass to 1 and tau-t to 2.0. ref-t was set to 120.27236 as suggested from the manual. However, the system was only stable when timestep was 0.0005, far smaller than the timestep used in the reference paper, 0.01. I used the same initial configuration and simulated with the same parameters in LAMMPS but the timestep could reach 0.006. I simulated the same system in HOOMD-Blue where the timestep can reach 0.01.

I’m struck by this result. I have carefully reviewed my GROMACS input files but have not identified any obvious errors. I would therefore like to know whether this instability of SD integrator is a known limitation of GROMACS’s algorithm for this type of system, or whether there is a likely error in my input file that I may have missed.

Please find attached my gromacs and lammps input files in the google drive link.

2

I see no reason why the SD integrator in GROMACS should be less stable than other implementations.

Have you checked that the energies at step 0 match between the different packages?

The google drive link requires access permission.

3

Dear Hess,

Thank you for your reply. I checked the initial energy in two packages, GROMACS gave a total potential energy of 7.27389e4, while LAMMPS gave 7.27376e4, which I think they are close enough.

I changed the access of the google drive link you should be able to view it now. Thank you!

I also attached the .mdp and .top here.

bd_eq_neu.mdp (2.4 KB)

bd100.top (7.7 KB)

4

I would think you should be able to use a time step of 0.001, but that depends on how well you initial configuration is equilibrated. But do I understand correctly that with LAMPPS the maximum time step is 0.006? In that case HOOMD seems to be the outlier.

5

Note that the friction is low, so the SD integrator should not significantly affect the stability. The maximum time step should be about the same with normal MD and a thermostat. That could help you identify the issue.