I don't know anything about details of your question, what Ermak and McCammon did etc. That being said "Is this right?" depends on what question you are asking. Say you want to solve: dx/dt = -lambda x w/ lambda > 0 so you discretize the left hand side: dx/dt -> (x(n+1)-x(n))/dt and the right hand side you replace with - lambda x(n+1) to find x(n+1) = x(n)/(1+lambda dt) = x(0)/(1+lambda dt)^(n+1) . Note that if we write t=n dt we have: x(t) = x(0)/(1+lambda t/n)^n which converges to exp(-lambda t) x(0) in the limit n-> infinity. But for finite t this correpsonds to the limit dt->0.

Now how large should dt be? If you only care about the long time solution you can take lambda dt -> infinity and you're system will converge to the correct long time answer. If you want to know the answer for finite times then you'll have to take a finite time step. If you want accurate solutions for intermediate time steps you'll have to have: lambda dt 2 the iterates diverge.

You're correct to question whether or not such a large time step is a good idea. Without more information, it's hard to say whether such an approach might make sense for a particular scenario. One nice reference that I think provides a good discussion of the trade-offs with coarse-grained models is http://dx.doi.org/10.1080/00268970310001640102.

To John: In case of stochastic differential equations the solution obtaied by a numerical method highly depends on the stpe size; relationships are much more compoicated than in the deterministic case.

Most deterministic problems are far more complicated than the one I discussed. To my knowledge the only reason to use an implicit schemes is because they side-step stability constraints imposed by explicit schemes. You will miss the fast dynamics if you do this. If your interest is the fast dynamics then obviously you'll get the "wrong answer". If your interest is the phase space density of trajectories over long times then you'll probably get something meaningful if the noise is handled properly.

My suggestion to your problem is a step adaptive method with an implicit solver.

I attach a paper with an approach that can be easily applied to any kind of algorithm.

In this way you start with a small time step and the adaptive method adjust it to have the required accuracy.

When you need to get a rapid change the step will be small, when the fast phase is finished the time step will increase to higher values. If you have high frequency driving force then the time step is limited by the frequency.

Hope it helps

@Pearson In their paper "Brownian dynamics with hydrodynamic interactions", Ermak and McCammon proposed a first-order integration algorithm for BD.

Actully, what I concern is the equilibrium mechanical property of crosslinked bead-spring chains so it may not involve dynamics.

Correct. If you can pose BD as an IVP, given that the solution will be stable and converges, an implicit time integrator can use larger time steps.