Hi!
I'm currently using the Langevin thermostat for a system of spherical particles, with the "omega" option enabled to thermostat the rotational DOF's. However, it seems to me that there's an erroneous factor in the random torque applied to the particles, as follows:
According to Stoke's law, the relation between the translational and rotational diffusion coefficients (D_t and D_r, respectively) for particles of diameter sigma, should be D_r = 3*D_t / sigma^2 in the overdamped limit. Moreover, the variance of the translational and rotational noise in Langevin / Brownian dynamics should be equal to 2*D_t and 2*D_r, respectively. In the "omega" part of fix_langevin.cpp, the equations are identical to those in the translational routine, with the substitutions m -> I (moment of inertia), f -> torque and v -> omega. Since for a sphere, I=(1/10)*m*sigma^2, this means that in the end D_r = 10*D_t / sigma^2, i.e. a factor of 10/3 larger than expected from Stoke's law. I have also confirmed this erroneous result from calculating the diffusion constants using the MSD and rotational correlation times, respectively, from a simulated dilute colloidal system.
I haven't looked into the "angmom" part of the code, but there might be an issue there as well. Could you please have a look at this and see if it's indeed a bug? Thank you!
Kind regards,
Joakim Stenhammar
Your analysis looks correct - I don't think we thought
about the connection to diffusivity when adding the
omega and angmom options. Let me check
with some colleagues. Reducing the applied torque by 3/10
would apply to only the random noise term, not the
drag term?
Thanks,
Steve
Hi again!
Sorry, but I didn't really think things through properly regarding the drag force - based on my discussions with Davide, I think this is how I think the code should be changed:
Rather than just making the transition m -> I in the case of rotational motion, care must also be taken to replace the friction coefficient gamma_t (equal to the inverse of t_period in Lammps) with its rotational counterpart gamma_r. Comparing with Stoke's law, one gets gamma_r = (10/3)*gamma_t, which in the end would lead to a factor of 3/10 inside the square root for the rotational noise factor, just as before, and will also change the drag force accordingly. I hope this reasoning sounds OK, and please get back to me with comments when you've looked at it! Thank you!
All the best,
Joakim
Hi again!
I just wanted to check whether you had time to look at this problem, and confirm the bug in the rotational thermostat? I also discussed a bit with a colleague who uses the thermostat of "fix rigid", and the rotational thermostatting there seems to have the same problem. Thank you for your help!
Kind regards,
Joakim
no - haven't gotten back to it yet - sorry
Steve
In the translational Langevin thermostat there
is a drag term with gamma1 that has m in it.
And a random term with gamma2 that has sqrt(m)
in it. We replaced m with moment-of-inertia
I for a sphere. You are saying that in both cases
I should be replaced with 3/10 I ? I.e. a reduction
of gamma1 and gamma2 by 3/10 and sqrt(3/10)
respectively?
Also, this does not affect the thermostatting properties
of Langevin on omega, only the variance (and thus diffusion)?
I.e. if you measured a rotational temperature from the
rotational energy (suitable normalized by 3 dof per particle)
we will get the same (correct) temperature with or without
this change?
I presume something similar is needed for angmom and
aspherical particles. Not clear if it would be a prefactor
like 3/10 on each component of the diagonal inertia
tensor ? That would at least reduce to the same effect
in the limit of aspherical -> spherical.
Steve
Hi Steve!
Sorry for being so slow, but I finally had time to have a look at the Langevin issue! It seems your analysis is correct, apart from that gamma1 and gamma2 should be increased rather than decreased. More explicitly, I've replaced the gamma1 and gamma2 definitions in the omega part of fix_langevin.cpp with
gamma1 = -(10.0/3.0)*inertiaone / t_period / ftm2v;
gamma2 = sqrt(inertiaone) * sqrt(80.0*boltz/t_period/dt/mvv2e) / ftm2v;
As you said, the change does not affect the rotational temperature, but only the rotational diffusion. I attach a plot of the orientational autocorrelation function for a system of spherical particles in 2D simulated using the new and the old codes; with the parameters used (t_period = 0.01, temp = 1.0, mass = 1.0, diameter = 1.0) the expected rotational diffusion constant is 0.03, which fits well with the data. The rotational temperature is reproduced correctly (fluctuating around 1.0) in both cases, although the timestep needs to be decreased slightly in the new code for the rotational motion not to blow up.
As you suggested, the corresponding change should be made in the "angmom" part of fix_langevin.cpp and in fix_rigid.cpp, for them to produce correct results in the spherical limit. As in the spherical case, I think it's best to just rescale gamma1 and gamma2 with a factor of 10/3 and sqrt(10/3), respectively, which will then affect each component of the inertia tensor. I haven't found any more routines dealing with thermostatting of rotational dof's, but if there are any, the same reasoning should be applicable there.
Finally, the fix you issued for the error in read_data with atom style dipole seems to have done the trick!
Does this sound reasonable to you? Please get back to me if you have more questions, or to confirm that things look OK!
Kind regards,
Joakim
chkdiff.pdf (32.3 KB)
1 Like
Hi Joakim - another user emailed about whether
it was good to add the 10/3 correction factors to the angmom
case, saying it will not be correct for non-spherical
particles. See details below.
Since the chief use of the angmom case is indeed
for non-spherical particles, e.g. ellipsoids, I’m thinking
he is correct. I.e. we should remove the 10/3 changes
we made to the angmom case and simply document
that the angmom case will not be correct for
rotational diffusion (but ok for thermostatting).
Comments?
Steve
Hi Steve!
Thanks for your email! This is correct - the changes to the Langevin thermostat case indeed does only apply to spherical particles. However, my idea was that it is better to use a prefactor so that the correct behavior is reproduced in the spherical limit, rather than one that is incorrect in all cases - so I don't see any reason why a prefactor of 1 is better than one of 10/3. How to choose the prefactor is not obvious, but either way I agree with Tomasz that it should be documented with a disclaimer that it only reproduces the correct rotational diffusion constant in the spherical limit (if the 10/3 factor is kept) or not at all (if it is changed back to 1). Someone else might have a better argument for keeping the original 1.0 factor, however!
Kind regards,
Joakim
I straddled the fence on this one. Fix langevin
now takes an angmom keyword in one of these 2 forms:
angmom no (the default)
angmom scale
where scale is a numeric value, e.g. 10/3 or 1.0 or whatever the user specifies.
So the user can use 10.0/3.0 if the particles are (nearly) spherical or 1.0
if not.
The doc page explains the details now. Will be in the next patch.
Thanks for the careful discussion/explanation. I wasn’t aware
that rotational diffusivity was sensitive to these pre-factors.
Steve