Dear Axel,
I have a problem about the Langevin dynamics in LAMMPS:
I'm simulating a system with Langevin equation with spatial dependent
damping coefficient gamma(x,y,z), so I cannot use the fix_langevin command
directly since gamma(x,y,z) is not a constant.
I try to solve the problem as follows:
1, use fix_addforce to add the spatial dependent damping force
-gamma(x,y,z)*v and the random force to the atoms;
2, use fix_nve to integrate the equation.
But I'm wondering whether this method will give the correct thermal
dynamics or not.
Do you know is it the correct way to do this? Thank you in advance!
Hi Ouyang, this looks correct formally, but your numerical accuracy may be lower because of two issues:
You are applying both Langevin forces together, and won’t be able to use integrators specifically optimized for the Langevin equation, like the GJF of fix langevin. But this is mostly an issue for “longer” time steps.
If gamma changes significantly over the spatial domain, it does make a difference whether you use gamma(x^{n}) or gamma(x^{n+1/2}) when updating the coordinates: x^{n} -> x^{n+1}. I’m not sure what is the best strategy for this, but for sure when using the standard integrator, you are assuming that gamma is step-wise constant instead of a continuous function. Again, any errors should go in proportion to the time step length.