This does appear to match up with what I would expect; the number of atoms is not a direct substitute for the degrees of freedom in the system (consider various rigid body restraints, shake etc). Should this not instead be using a value derived with “temperature->dof”, rather than the number of atoms?
Aidan can probably answer better. The virial
portion of the pressure is initially computed for
the entire system and the constraints are subtracted
out later (for the virial part, kinetic contribution
is different). Also, the barostat mass is an artificial
construct anyway, so if it is 10% different, it is
I have another question regarding the Nose-Hoover barostat implementation. Looking at the integrator split in initial_integrate() and final_integrate(), I’m a bit confused.
It looks like the code considers the “symmetry” in the order of function calls to be equivalent to the “symmetry” which appears when splitting short-time approximations to Liouville propogators via Trotter’s approach.
For example, the nh_v_press() routine scales the particle velocities using the previously calculated barostat velocities, including the MTK energy correction term. The first time this routine is called ( in initial_integrate() ), the internal barostat variables have just been updated. However, the second time this routine is called ( in final_integrate() ), it is still using the “old” barostat variables.
If this routine is genuinely independent of the intervening changes to particle position and momenta, would it not make more sense to combine these two calls into a single application of the scaling factor? If the routine is dependent on changes to the particle momenta and positions (perhaps indirectly, via their effects on the barostat variables due to changes in instantaneous pressure components), then this looks very suspicious.
Likewise, the two applications of the remap() routine appear redundant: this routine rescales the simulation cell according to the barostat velocity, which does not change between calls to remap(). Calling it twice looks like it might be a NOOP, other than the fact that it evolves the barostat position (which is only really useful for checking for a conserved extended system).
Can anyone help me to understand the design of this part of the code a little better?
This formula for barostat mass follows the expression for the time constant commonly used in the literature. The objective is to translate the specified value of p_freq into a barostat trajectory that exhibits roughly the same time constant. In this context, differences on the order of O(1/N) are of no consequence.
Just because you are confused and/or feeling suspicious does not necessarily mean the code is wrong, or even poorly designed. If you look at the definition of FixNH::remap() you will quickly see that it is never a NOOP. The symmetry of the function calls is not just for aesthetics, although that would be a good enough reason by itself. The second call to nh_v_press is using the same omega_dot values as the first call, but the velocities have been incremented twice in between.
Just because you are confused and/or feeling suspicious does not necessarily mean the code is wrong, or even poorly designed
Very true - and that’s why I didn’t say the code was wrong and/or poorly designed, I was just trying to get a better handle on it!
If you look at the definition of FixNH::remap() you will quickly see that it is never a NOOP.
NOOP was a poor choice of wording on my part, you’re right. I’ll try again:
FixNH::remap() updates the simulation cell using the current cell momenta, and it’s called immediately before and after the nve_x() method. The nve_x() method does not alter the simulation cell, or the cell momenta: hence, I was wondering why this update occurred in two calls to remap() rather than a single call which did the full update (with suitable modification of the “dto” variables).
Sorry for the brusqueness, I just wanted to make sure you were actually looking at the code before shooting off questions. To answer your restated question, combining the two calls to remap() into one would not give exactly the same trajectory, for the same reason as mentioned previously i.e. the positions are incremented in between the two calls. This is not a big effect, but it is better to remain consistent with the original equations.