regarding nphug implementation

Dear all,

Dear all,

I am using “fix nphug” to study shock wave energy dissipation. My command reads as below:

fix myhug all nphug temp 1.0 1.0 10.0 x 78953.861 78953.861 1000.0 (units real)

When I output “f_myhug” using the thermo_style command, does this “f_myhug” correspond to “E(H) = E0 + 0.5(P+P0)(V-V0)” or “E(H) - E”?

As f_myhug is zero for the first step, I am assuming it is “E(H) -E” as E0 = E for the first step and V = V0.

Can anyone please tell me if it is implemented like this?

From the original publication, I gather that the parameters used for controlling the ergostat and the barostat play a very important role in these simulations. From a previous post on the forum, I gather that these variables are linked to the “tdamp” and “pdamp” parameters in lammps implementation.

Can anyone also tell me if in nphug, the conserved quantity is “E(H)” or “E(H) - E” with time?

I would like to see if I am conserving the right quantity with my chosen set of tdamp and pdamp parameters.


Karthik Guda

Purdue University

The documentation is fairly clear about all this:

"These fixes compute a global scalar and a global vector of quantities, which can be accessed by various output commands. The scalar value calculated by these fixes is “extensive”; the vector values are “intensive”.

The scalar is the cumulative energy change due to the fix.

The vector stores three quantities unique to this fix (Delta, Us, and up), followed by all the internal Nose/Hoover thermostat and barostat variables defined for fix npt. Delta is the deviation of the temperature from the target temperature, given by the above equation. Us and up are the shock and particle velocity corresponding to a steady shock calculated from the RH conditions. They have units of distance/time."

The quantity you are looking for is related to f_myhug[1], which is defined by the equation given. The scalar quantity f_myhug is something else entirely. Just as for the various Nose-Hoover barostats and thermostats, adding this quantity to Etot = PE+KE, gives a quantity that is approximately conserved by the time integration scheme, depending on how small a timestep is used.


Dear Aidan Thomson,

Thank you for your reply.

Can you please tell me if conserving the quantity (f_myhug + Etot (KE+PE)) with time equivalent to achieving the “Hugoniostat energy condition (E(H) - E) ~ 0” for a steady state shock simulation? Or should I monitor the quatity “f_myhug[1] (Delta)*N(dof)*kb”, based on the equation given as a function of time in order to see if the hugoniostat energy condition is met? I also see that in addition to the timestep, changing “tdamp” parameter results in different behavior for the quantity (f_myhug + Etot) with respect to time. So, can I assume that both of them can be varied so as to conserve f_myhug+Etot?



There are two separate issues here. First, one would like to verify that the time discretization of the equations of motion is sufficiently accurate. Monitoring conservation of f_myhug + Etot (KE+PE) is one way to do that. Second, one would like to know whether the system has evolved to a stationary state that is sampling a distribution of positions and velocities that are consistent with the Hugoniot jump conditions. Monitoring f_myhug[1] is one way to do that. For the most part, these two issues are independent or each other.

Dear Aidan Thompson,

Thank you for your reply.