For #1, you're thinking about it the wrong way. Usually, one can reasonably assume that for small enough dt, the error rate on any observable should asymptotically approach a constant. What this suggests is that you should run your simulation at dt_1 and show that your observable is insensitive to changes in dt by running it at dt_2 < dt_1. In other words, don't try to calibrate dt by waiting till you get the expected result - for a given dt that _gives the correct result (within what you define to be an acceptable error) - show that the answer you get is insensitive to perturbations in dt.

#2 Without knowing the underlying code, it's hard to say. However it does look like the system has not fully equilibrated yet. Any thermodynamic equilibrium values you compute at this point would be incorrect if that is the case.

#3 If the timestep does not give you a known expected benchmark value - you are working on a different system. Remember, numerical integration is only an approximation to the underlying diff eq's. If a known result is different, how could you possibly expect an unknown result to be the correct value?

Ok thank you very match

You should first estimate the smallest time scale in the problem (inverse friction, bond oscillation frequency) and set your time step such that dt < smallest time scale. For bead-spring chains the relaxation time grows like N^(2 nu + 1) with Rouse dynamics, so it takes a very long time to equilibrate long chains N >> 1. All these technical aspects can be found in the many, many papers published on this topic and also standard books such as Allen & Tildesley.

@Travis Thank you for tellin me these. I think I have better understanding of this now. I used smaller time step and got consistent results.

@Tapio Thank you. yes now I realize it's because of the long relaxation time of Rouse chain model. And I have another puzzle: how to estimate the bond oscillation frequency in brownian dynamics neglecting the inertia term?

You're welcome! If you don't have the inertial term, the bond oscillations are overdamped. The microscopic time scale is then inverse friction.

@Tapio So do you mean the stiffness of the spring doesn't have effect on determining the time step? this seems not to be the case as seen in my simulation.

Yes, of course the spring constant affects your time step - the larger the forces i.e. the steeper your V(x-x'), the smaller the time step has to be to prevent unphysical overlaps. There's a huge literature on this - one can also use adaptive (smaller) time steps when overlaps occur. Search Google!

@Tapio OK, thank you again!

What numerical integrator do you use?

If you pick a time step dt that is too long, I would expect the variance of thermal fluctuations to increase, not necessarily the mean. What you plot looks to me very much like a system that equilibrates fast (within t~2000) but that does not conserve total energy over time.

A simple Euler integrator, for example, will constantly overestimate positional displacement, because it only estimates repulsive forces at the far end of the interaction, thus leading to a too strong particle overlap. When the particles separate, it will equally overestimate the particle forces, thereby increasing the kinetic energy. When integrating the equations of motion, what you will like to use is a symplectic integrator that keeps the phase space volume (total system energy) constant. An easy to implement symplectic integrator is the Verlet integrator.

Of course, if you are already using a suitable integrator, just ignore my input...

It is possible to use a simple step adaptive criterion that can be easily adapted to any kind of algorithm. I attach the paper. The method is based on guessing the truncation error comparing two different solutions. Obviously you need the proper algorithm.

By the way you can use the energy as the quantity used to check the time step to limit the energy growth. A very simple symplectic algorithm is the lipfrog method, which solve alternatively the equation for velocity and that for position.