Today a colleague and I noticed something that has got us stumped. When running a simulation with a simple polymer chain of point particles with harmonic bonds, the observed average bond length is systematically larger than the expected equilibrium bond distance.
Attached are a minimal input script that produces this problem along with a data file it includes. There is no pair-pair interaction, just bonds. Expected behaviour is that the average bond distance should fluctuate around the equilibrium distance specified in the input script (1.0 in this case), but it fluctuates about a slightly higher value (1.019). Adding more chains reduces the fluctuations about this value but does not change the average. Many test cases we tried (more/fewer bonds, as low as one bond per "poly"mer, many polymers, other equilibrium distance, other spring constant) all reproduced this problem.
I am currently using the LAMMPS version of August the 29th, 2014.
I am really not sure why this happens, is there something we are missing?
We think we have discovered the origin of this observation today, and it turns out that it is not actually a problem but rather physics.
Steve, your assumption that a harmonic bond should not have an average of 1.0 exactly for finite temperatures is correct, and yes, the stretch increases with temperature.
The deviation occurs because for each indentation of the spring there is an extended state with more possible configurations, corresponding to all points on sphere surfaces of some radius r_0 +/- dr, with r_0 the equilibrium bond length and dr the indentation/extension. Both the extension and indentation have equal energies of k * dr^2, but there are slightly more states for the r_0 + dr configuration.
Simply plugging this into a Boltzmann-factor and doing all the statistics already shows that the expectation value for the actual bond length is larger than r_0. The expression is a bit complicated, for large enough kappa / kT it reduces roughly to = l_0 + kT / (k l_0), with l_0 the equilibrium bond length and k the spring constant (as given in the input script), and kT the thermal energy. The expression we get agrees at least in the first decimal with the deviation from l_0 we observe in the simulations, so this seems like the proper explanation of the observation.