Hi Espressies,
Simulating a single polymer chain of the Kremer-Grest model, we observed that the centre of mass "diffuses" in narrow a slab, with a normal vector approx. (1,1,1). In order to see that, one needs to run about 10^8 LJ time units (~10^10 time steps, ~1cpu-day), plot the centre of mass trajectory and tilt it in an appropriate angle. With shorter trajectories it's not obvious if the bias is in the trajectory or in the observer's eye.
I could partly fix the problem replacing the call to l_random() (ran1 RNG) in d_random() of random.hpp with the built-in rand() function. The trajectory might still be biased, but at least it's not a slab.
Attached files: sample tcl script, gnuplot script to plot what I have described and its source data: centre-of-mass and monomer msds (average over 20 runs), com trajectory with the l_random() and with rand() RNG.
I think that we should replace the l_random() by some other RNG. Surely, rand() is not a good choice either, it was just the first hack I tried. Any suggestions / opinions are welcome. I intend to investigate it further: test other RNGs and check if/how they affect dynamical and thermodynamic properties.
Some more detail of what I learned:
I asked a student to simulate the Kremer-Grest polymer in Langevin thermostat and compute msd(t), from 20 independent starting configurations, to get 20 independent msds's and estimate the error. Plotting msd(t)/t, it first converges to approx. 2/N, as it should for Rouse-like dynamics, but at t~3*10^4, it drops to ~3/4 of this value, while the error is much smaller. We suspected a mistaken input or a subtle correlator bug, so I repeated the simulations with my own script, various chain lengths, box sizes, gammas, timesteps, correlator parameters, initial conditions, and the result was robust. Finally, I printed the centre of mass coordinates and asked a colleague to compute the msd(t) using his own independent tool. Plotting the output trajectory revealed the problem. The tilt of the plane appears robust with respect to initial conditions, random seed and chain length. I have not yet checked other parameters, since other runs did not write the trajectory. I also checked the original Grest+Kremer paper (Phys.Rev.A,33,3628,1986), and their msds also change slope on long times (g3 in Fig.3), but they have not specified their RNG.
Looking forward to your comments.
msds.tar.gz