Hi everyone,
in my bulk system I observe huge pressure oscillations. It is my understanding that these are to be expected for such a system and ought to be damped by the drag parameter.
Unfortunately, I failed to understand what exactly the drag parameter does and how it affects the NH barostat. From older mailing-list replies and the manual, we can learn that drag is an ad-hoc fix to damp oscillations in the NH barostat, which alters the NH formalism. Unfortunately, it remains unclear how exactly the formalism is changed.
-
I did not find any further documentation or a reference that would explain what is done and how this affects the results in detail, did I miss it?
-
In case there is no such thing in writing, I guess the answer has to be looked for in the code. Looking through the source files, I am guessing that this “dragging” is done in fix_nh.cpp. A pdrag_factor is created that seems to remain 1 without drag (i.e.,drag=0.0). Otherwhise, 1 is reduced by a factor depending on the time step, Pdamp?, drag, and the number of chains (e.g., line 689).
Next, pdrag_factor is multiplied to etap_dot and omega_dot, which should be the barostat velocities (omega_dot) and the thermostat velocities belonging to the barostat (etap_dot).
From this, I assume that the pdrag_factor reduces the resulting change that the heat baths would normally (i.e., without drag) induce, by scaling down the velocity changes. Could anyone confirm this?
Next, how can we interpret this? I guess it is directly evident that one may loose energy conservation by this procedure. As a simple example: For the next step, it is determined that the bath is supposed to give energy to the system, by increasing the velocities by delta v, however, due to drag, the systems velocities only obtain a change of delta v’, where delta v’ < delta v. For energy conservation, we would have to add up E(System) + E(Bath) + a factor that is related to the difference delta v’ - delta v. If LAMMPS is not accounting for this additional factor, energy conservation is lost. What does this mean for our results?
Dr. Julian Gebhardt
Materials Modeling
Business Unit Materials Design
Fraunhofer Institute for Mechanics of Materials IWM
Woehlerstrasse 11 | 79108 Freiburg | Germany
Phone + 49 761 5142-255 | Fax + 49 761 5142-510
julian.gebhardt@…7756… | www.iwm.fraunhofer.com
Your point about drag is mentioned in the documentation of fix_nh:
“In some cases (e.g. for solids) the pressure (volume) and/or temperature of the system can oscillate undesirably when a Nose/Hoover barostat and thermostat is applied. The optional drag keyword will damp these oscillations, although it alters the Nose/Hoover equations. A value of 0.0 (no drag) leaves the Nose/Hoover formalism unchanged. A non-zero value adds a drag term; the larger the value specified, the greater the damping effect. Performing a short run and monitoring the pressure and temperature is the best way to determine if the drag term is working. Typically a value between 0.2 to 2.0 is sufficient to damp oscillations after a few periods. Note that use of the drag keyword will interfere with energy conservation and will also change the distribution of positions and velocities so that they do not correspond to the nominal NVT, NPT, or NPH ensembles.”
Dear Stefan,
thank you, I was, of course, aware of that section in the manual. I am wondering about the specifics:
“drag keyword will damp these oscillations, although it alters the Nose/Hoover equations” - how does it change the NH equations?
“A non-zero value adds a drag term;” - how does the drag term look like?
“Note that use of the drag keyword will interfere with energy conservation and will also change the distribution of positions and velocities so that they do not correspond to the nominal NVT, NPT, or NPH ensembles.” - what does this mean for results obtained with a drag term?
Best,
Dr. Julian Gebhardt
Materials Modeling
Business Unit Materials Design
Fraunhofer Institute for Mechanics of Materials IWM
Woehlerstrasse 11 | 79108 Freiburg | Germany
Phone + 49 761 5142-255 | Fax + 49 761 5142-510
julian.gebhardt@…7756… | www.iwm.fraunhofer.com
For the shape of the drag and stuff, the answer is in the code. I am not familiar with fix np* so I cannot really help there. I do agree that the documentation could have been a little more explicit there.
As for the sampling, it means exactly what it says. Results obtained with a drag term will not sample the “correct” ensemble, i.e., they will not correspond to the true np* ensemble. It is hard to say how they deviate because that will depend strongly on the type of system you simulate. Maybe in your case it won’t be so bad. The only way to know for sure is to simulate with and without drag, perform sampling for both, and seeing if the parameters you are interested in are significantly different.
When dealing with constant pressure simulations, it might be helpful to use drag (or even a way too dissipative barostat like Berendsen) to settle your system into the right average pressure and then sampling with a proper barostat without drag once you reach a steady state. That in itself might already significantly reduce the fluctuations. If they are still too bad, you can also experiment with the number of barostats through the pchain option. As far as I know, barostat chains will sample the right ensemble.