Dear Lammps users,
I am running some simulations in Lammps using Nose-Hoover temperature thermostat. I have a question about how Nose-Hoover dynamics is implemented in the Velocity-Verlet algorithm. As I have found in the literature that in the Nose-Hoover dynamics, acceleration depends on velocity (which is not the case for Velocity-Verlet algorithm). Therefore, the new velocity appears on both side of the equation while solving it for the velocity. My question is how does Lammps solve the velocity and acceleration when Nose-Hoover thermostat is selected. Does Lammps solve it by an iterative scheme like predictor-corrector algorithm, instead of Velocity-Verlet algorithm?
I would greatly appreciate any suggestion or reference.
Thanks,
Ishraq Shabib.
In the Nose-Hoover (NH) scheme, the velocity V of a given particle obeys the equation
dV/dt = F / m - eta * V,
where F is the total force on that particle (due to other particle interactions and/or external forces) and eta is a variable representing the interaction of the system with a heat bath.
Integrating equations of motion involves discretizing time. This gives for the equation above something like
V(t+dt) = V(t) + [F(t) / m - eta(t) * V(t)] * dt.
Hence the velocity at time t+dt is dependent upon the velocity at time t. (the “new velocity” only appears on the left hand side of the equation)
The NH integration scheme is actually more complicated than the second equation I wrote down, as one needs to worry about time-reversibility and accuracy of the integrator in orders of dt. If you are interested in these details I would take a look at “Understanding Molecular Simulation” by Frenkel and Smit. There you will find discussion on how integration schemes are performed in general, as well as a few sections specifically devoted to NH dynamics.
Hope this helps.
Andy
Ishraq,
No, an iterative scheme is not required.
The procedure LAMMPS uses looks like velocity verlet (½ v update, full x update, f update, ½ v update).
Take a look at the source code of fix nvt to see the details.
See also:
http://lammps.sandia.gov/doc/fix_nvt.html
Hoover, Phys Rev A, 31, 1695 (1985).
Paul
Dear Paul and Andy,
thanks for your helpful comments. I have read the book “Understanding Molecular Simulation” by Frenkel and Smit for numerical implementation of Nose-Hoover dynamics in velocity verlet algorithm. The algorithm looks like what Paul suggested -
(½ v update, full x update, f update, ½ v update).
However, I’m still confused about the way the last ‘1/2 v update’ is done. according to the book by Frenkel (page 387 and 391)
v(t+dt) = v(t+dt/2) + [ f(t+dt)/m - eta(t+dt) v(t+dt) ] * dt/2 ---------- (1)
eta(t+dt)= eta(t+dt/2) + [sum (m* v(t+dt)^2) -gT ] dt/2Q ------------(2)
In both of these two equations we see the “new velocity v(t+dt)” term on the right hand side of the equation, and the book mentioned Newton-Raphson scheme to solve them numerically.
I then looked at the source code and the documentation page for fix_nvt. I found that the way the thermostat works is-
“the current temperature is calculated taking the bias into account, bias is removed from each atom, thermostatting is performed on the remaining thermal degrees of freedom, and the bias is added back in.”
My understanding is when BIAS is removed from each atom the advection velocity is removed, but the thermal velocity remains. Does that mean when BIAS is removed it somehow influences the v(t+dt) term on the right hand side of equation (1) & (2) above?
or as Andy said "… the velocity at time t+dt is dependent upon the velocity at time t. (the “new velocity” only appears on the left hand side of the equation (1). "
I’m sorry for this long email. I would really appreciate any suggestion to clarify my understanding.
Thanks,
Ishraq Shabib.
Ishraq,
You have raised two separate questions
- Does the LAMMPS NVT thermostat update velocities appropriately?
- Does the LAMMPS NVT thermostat handle net motion appropriately?
The answer in both cases is yes.
It is true that some codes treat the second Verlet update as an implicit equation that must be solved iteratively. However, it seems that this is overkill and also breaks the nice time-symmetry of the standard velocity-Verlet scheme. We have followed the work of Martyna, Tuckerman and Klein, who have shown that explicit equations, when solved in the correct order, sample the canonical ensemble, which is the whole point of NVT.
The second issue of net motion (bias) is necessarily more ad hoc, but the scheme that is in LAMMPS is pretty reasonable.
If you have a specific problem let us know. If you are just curious about the details, the best thing to do is read the ultimate documentation, the code itself.
Aidan