Dear all,
I have just finished my Bachelor’s and I am in a research internship on ABP. We are currently reviewing the dumbbell model with bond_style FENE (bond_style fene command — LAMMPS documentation) , in particular, studying the diffusion coefficients for different sizes of the dumbbell. The relation between these coefficients and this parameter has already been derived theoretically (see; Phys. Rev. E 91, 062124 (2015) - Rotational and translational diffusion in an interacting active dumbbell system, eq 30 and 43), but we want to check numerically that this is still true in the heavily over-damped limit.
So far, we had done our simulations with fix nve+fix langevin
fix 1 all nve
fix 2 all langevin $T T (1/v_gamma) 457356
and we choose arbitrary \sigma, R and k (see: ) such that U'(a)=0, where a is the distance between the two particles and U is the sum of FENE and Lennard Jones. This allows us to determine an \epsilon=\epsilon(\sigma, R, k,a). With the parametrization, we achieved a reasonable distribution of the a parameter during the simulation and matched theoretical predictions for the coefficients.
However, once we simulate with fix Brownian:
fix 1 all brownian $T 12908411 gamma_t 10
we do not obtain an acceptable distribution for the a parameter: it tends to be largely asymmetric and towards smaller a values than set. We have tried many combinations of R, k and \sigma but we are not satisfied with our results. We have even changed the parametrization to choose arbitrarily \epsilon instead of \sigma and find \sigma=\sigma(\epsilon, R, k, a) such that U'(a)=0 but this has not provided satisfactory distributions for a parameter. Also, we usually obtain the warning for too long bond, but we can control this with R. But some values of R are acceptable for fix Langevin and not Brownian.
For the moment we simulate with no interaction between dumbbells. a ranges from [0.1,1] and k=0.1, R=0.55, \sigma=0.486 (computed with the last method explained) and a=0.5, \gamma_t=10, T=0.1 and \Delta T=10^{-6} we have obtained the best results, but I wonder
- is there any fundamental difference between the two fixes that makes the bond_style act so differently (obviously they are different dynamics and regimes)?
and
2)how we can choose the parameters such that we obtain narrow and symmetric distributions of the parameter a.
I hope I have provided enough detail.
Thank you all very much,
M.