Hello,
I believe I found an error in the Bussi thermostat implementation in LAMMPS. (Actually, I think the error is in the equation of motion given in the Bussi paper). It’s not very serious, but for a low number of degrees of freedom some errors could occur.
It was easier to type it in Latex than in email, so I’m attaching it here. I can make the fix and push it to the repo, but I figured I’d check first that the fix is wanted and that I’m not making an error myself.
Hello,
I believe I found an error in the Bussi thermostat implementation in LAMMPS.
(Actually, I think the error is in the equation of motion given in the Bussi
paper). It's not very serious, but for a low number of degrees of freedom
some errors could occur.
please note that LAMMPS contains two (2!) thermostats following papers
of giovanni bussi, fix temp/csvr and fix temp/csld.
both have been implemented following the reference code provided on
giovanni's home page. fix temp/csvr also contains fragments provided
by paolo raiteri. there should be some reference in the mailing list
archives.
It was easier to type it in Latex than in email, so I'm attaching it here. I
can make the fix and push it to the repo, but I figured I'd check first that
the fix is wanted and that I'm not making an error myself.
you probably should first and foremost discuss this with giovanni
bussi directly.
we're happy to include code improvements. for obvious bugfixes, or new
features, that is straightforward. if you are changing the
implementation, some kind of external confirmation that your change is
useful and meaningful is usually much appreciated.
axel.
Good idea to contact Dr. Bussi. I was looking at his webpage to find his contact info, and on it he actually mentions this:
“(note) There is a typo in the Bussi, Donadio and Parrinello JCP (2007) paper in the instructions for calculating the sum of Gaussian numbers using a gamma distribution - many thanks to Francois Gygi for pointing it out! If you follow exactly the instructions in the paper you will end up in a wrong implementation. On the other hand, the routines that you can download from this website are working properly. In case you need help, feel free to contact me.”
So multiplying by 2 is what Bussi did himself to fix the issue.
I actually can’t find a better implementation of a gamma distribution with beta=0.5 than one in which one takes the gamma distribution with beta=1 and multiplies it by 2. Numerical Recipes: The Art of Scientific Computing only gives a recipe for beta=1, and per https://en.wikipedia.org/wiki/Gamma_distribution#Generating_gamma-distributed_random_variables, there might not be a better method for low Nf values.
So best to leave this be.