# [lammps-users] How many timesteps are enough for equilibration? Checking the sampling quality on-the-fly

Hello Lammps community,

First of all, I would like to apologize if my question is not directly related to Lammps.

I have a possibly simple question about the sampling quality, not the Lammps itself. To equilibrate my system (A bead-on-spring linear polymer of length 2000 monomers crowded by WCA beads in a cylindrical condiment), I use Langevin thermostat and fix NVE/limit with dt=0.01 and 10e6 timesteps.
I then run the system with dt=0.005 for 5e7 timesteps in the sampling phase. I dump the configurations of the system every 1000 timesteps in the sampling phase.

I use the block-averaging to estimate the statistical uncertainty in my system, as described the standard textbooks on computer simulations (Allen and Tildeskey, Rapaport, Frenkel & Smith). If I use the dumped configuration to estimate the statistical uncertainty, I have 5e4 configurations instead of the original 5e7 ones. Because of this reduced number of configurations, the number of block operations that can be done on 5e4 configurations is significantly smaller than the number of block operations can be done on the original 5e7 configurations. Does this cause any problems in reaching the plateau that should be seen in block-averaging method?

Let’s give an example, I use the dumped 5e4 configuration to compute the radius of gyration of the polymer and then use the block-averring method to find the statistical uncertainty in the mean radius of gyration of the polymer. According to the attached plot, the system is not run enough to reach equilibrium since there is no plateau in the statistical uncertainty vs the number of block operations. Has this failure in reaching the equilibrium happened due to the reduced number of configurations (5e4 instead of 5e7)? Or, even if I used the total number of configurations, I would not still see the plateau? If I used the all the configuration, I had 25 block operations instead of 14 operations - the vertical line for each data point shows the standard error for that data point.

According to the textbooks and the original paper (https://aip.scitation.org/doi/abs/10.1063/1.457480?casa_token=9nGutL6jMI8AAAAA:zr2RSixcQCFqa4_40IUvBZsOB_kAoprc8gvgi9VJPPlgUQXcmVXzFMgUPgyUtrINyusmkw5G3Oo) on the block-averaging technique, one can see the plateau after 6-8 block-operation.

This experience triggered a more fundamental question for me: How can I check whether the system has reached equilibrium or not? Can I do this check on-the-fly in lammps?

Since my question is not about a lammps problem, I do not include a log file.

Kind regards,
Amir

[…]

This experience triggered a more fundamental question for me: How can I check whether the system has reached equilibrium or not? Can I do this check on-the-fly in lammps?

the answer is - as with many of these kind of questions - it depends. there are no guarantees. comparing blockwise averaged results with bisectioning is a good approach, but it cannot tell you about “rare events” or how shallow or deep your potential hypersurface is and whether there are other, relevant minima in phase space that you don’t get to sample because it requires a very large energy fluctuation to get there. of course, this strongly depends on the nature of your system.

to give an example: both the CHARMM and GROMOS force fields are tested to represent available crystal structures of proteins well. but they differ in some design choices, for example in how “stiff” dihedral potentials are. CHARMM is much more stiff than GROMOS. That means with CHARMM any MD simulations “look” converged more easily and with GROMOS it takes more simulation time to get results that “look” converged. but then again CHARMM structures are more often “meta stable” and require more time or extra simulation tricks to boost sampling the “dihedral phase space”.

in the end, except for trivial systems, you will never know for certain. sometimes experience helps, but even the most experienced people can overlook something.

axel.