Hi, all-
I’m using MD to simulate FCC Lennard-Jones. I am not getting consistent values between compressibility determined by finite difference and by fluctuation formulas. I tried it in reduced units, was getting inconsistent results, and switched to real units to make it easier to compare numbers. I’ve also compared to GROMACS, where I get much better consistency. This appears to be an underlying issue with LAMMPS not getting the proper Boltzmann distribution for the NPT ensemble. I wonder if there is some subtle setting I’m getting wrong?
Some info:
Formula for finite difference used = - <Vhigh>-<Vhlow>/(((Phigh-Plow) (0.5(<Vhigh> + vlow>)))
Formula for finite difference by fluctuation: std(V)^2/kBT*<V>, kB = 0.1380 bar/nm^3 / K.
LAMMPS:
T P <V (A^3)> dV/<V>
44 797 36833.043 0.002410
44 897 36751.109 0.002408
compressibility from finite difference = 2.2565e-05 atm^-1
compressibility from fluctuation (at 797 atm) = 3.5237e-05 atm^-1
compressibility from fluctuation (at 897 atm) = 3.5108e-05 atm^-1
compressibility from fluctuation (average) = 3.5173e-05 atm^-1
GROMACS:
T P <V (nm^3)> dV/<V>
44 797 37.33091 0.002465
44 897 37.19381 0.002427
compressibility from finite difference = 3.6792e-05 bar^-1
compressibility from fluctuation (at 797 bar) = 3.7367e-05 bar^-1
compressibility from fluctuation (at 897 bar) = 3.6073e-05 bar^-1
compressibility from fluctuation (average) = 3.6720e-05 bar^-1
I expect the fantastic agreement of GROMACS to be somewhat a statistical fluke, but it is consistent.
I’ve also run ‘checkensemble’ (https://github.com/shirtsgroup/checkensemble) on the volume fluctuations. In this case, I got for the maximum likelihood analysis of the slope of log P_1(V)/P_2(V):
GROMACS:
LJ_44_897.mdp (1017 Bytes)
LJ_44_797.mdp (993 Bytes)
in.real.44_797.LJ (682 Bytes)
in.real.44_897.LJ (683 Bytes)
Aidan (CCd) can likely give some insight on this.
One note: in your scripts, you are using a time constant
of 5.0 with real units on both the thermo and barostats.
That’s 5 fs in real units, which is quite short and could
affect the fluctuations. More typical would be 100 fs
for the thermostat and 1000 fs for the barostat.
5.0 in LJ units would be fine, since it is ~1000 timesteps.
Steve
Sorry, I didn’t see that your timestep is 0.005
in the scripts. I assumed it was 1 fs (the default).
That’s a 0.005 fs timestep for liquid Ar or the like?
Steve
Following up on Steve Plimpton’s pointing out that I was assuming ps as the real time unit, rather than fs, I tried to rerun to make the simulations more comparable. However this made no difference: it was statistically identical, though it’s numerically even clearer that there is an issue, since the number of uncorrelated samples increased with the 1000x longer timestep. I changed the time for temperature and pressure pcoupling to 400 and 2000 time steps to match the GROMACS simulations more precisely. I’m attaching the new files.
Note that I independently, added ‘mtk yes’ to the previous simulations (i.e. changed “fix id1 all npt iso 797 797 5.0 temp 44 44 5.0”, to “fix id1 all npt mtk yes iso 797 797 5.0 temp 44 44 5.0”), which resulted in literally no change – it was identically numerically out to at least 50K steps.
Also, I verified that the average pressure estimator output by LAMMPS is statistically correct (for 797 simulation, pressure is 796.6 +/- 0.9 . So the average pressure is right (I didn’t expect that to be wrong, so good 
.
Data below:
BEFORE (time step really small):
T P <V (A^3)> dV/<V>
44 797 36833.043 0.002410
44 897 36751.109 0.002408
compressibility from finite difference = 2.2565e-05 atm^-1
compressibility from fluctuation (at 797 atm) = 3.5237e-05 atm^-1
compressibility from fluctuation (at 897 atm) = 3.5108e-05 atm^-1
compressibility from fluctuation (average) = 3.5173e-05 atm^-1
Validation of NPT ensemble looking at slope of log P1(V)/P2(V)
in.real.44_897.LJ (689 Bytes)
in.real.44_797.LJ (688 Bytes)
I discovered an inconsistency in the cutoffs – hold off on thinking about this until I’ve had a chance to regenerate the data!
Best,
“mtk yes” is the default for that fix npt option.
So that’s why specifying it explicitly didn’t
change anything.
Steve
This is not a LAMMPS question, but an NPT-dynamics question. If you read the literature on this subject, you will see that volume fluctuations from any Nose-Hoover-style NPT simulation do not reliably converge to the true isothermal compressibility. This is particularly true for crystals. The reasons are unclear, but it probably has to do with effective non-ergodicity, in which the system gets trapped on very long-lived manifolds in the space of volume fluctuations. If GROMACS has found some magic fix (it’s hard to tell from the thread if this is was observed or not), I would love to know what it is.
For a very thorough study of this problem from 33 years ago, see Sprik, Impey, and Klein, PRB 29 4368 (1984).
Aidan
PhysRevB.29.4368.pdf (503 KB)
So, It turns out that as I converted over the files from reduced to real units, I had made an error and did not have matching cutoffs. It turns out that I was able to get both real units and reduced units working. In the case of reduced units, with T = 0.3, and P = 2.0 and P = 2.3. I get the following result for the compressibility calculations for LAMMPS:
compressibility from finite difference = 1.4006e-02
compressibility from fluctuation (low) = 1.4084e-02
compressibility from fluctuation (high) = 1.3869e-02
compressibility from fluctuation (average) = 1.3977e-02
And for the checkensemble check of volume distribution slope: