This is not a LAMMPS question, but an NPTdynamics question. If you read the literature on this subject, you will see that volume fluctuations from any NoseHooverstyle 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 nonergodicity, in which the system gets trapped on very longlived 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
Sandia National Laboratories
PO Box 5800, MS 1322 Phone: 5058449702
Albuquerque, NM 87185 Fax : 5058457442
Email:athomps@...3... Cell : 5052181011
From:
Steve Plimpton <sjplimp@...24...>
Date: Friday, July 21, 2017 at 3:45 PM
To: Michael R Shirts <Michael.Shirts@...780...>
Cc: lammps/lammps LAMMPS Users List <lammpsusers@lists.sourceforge.net>
Subject: [EXTERNAL] Re: [lammpsusers] Issue with pressure control not giving the correct fluctuation properties with LennardJones solids
"mtk yes" is the default for that fix npt option.
So that's why specifying it explicitly didn't
On Thu, Jul 20, 2017 at 3:52 PM, Michael R Shirts <Michael.Shirts@...780...> wrote:
I discovered an inconsistency in the cutoffs – hold off on thinking about this until I’ve had a chance to regenerate the data!
Best,
~~~~~~~~~~~~~~~~
Michael Shirts
Associate Professor
michael.shirts@...780...
http://www.colorado.edu/lab/shirtsgroup/
Phone: (303) 7357860
Office: JSCBB C123
Department of Chemical and Biological Engineering
University of Colorado Boulder
On 7/20/17, 1:54 PM, "Michael R Shirts" <Michael.Shirts@...1657...> wrote:
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.2565e05 atm^1
compressibility from fluctuation (at 797 atm) = 3.5237e05 atm^1
compressibility from fluctuation (at 897 atm) = 3.5108e05 atm^1
compressibility from fluctuation (average) = 3.5173e05 atm^1
Validation of NPT ensemble looking at slope of log P1(V)/P2(V)

Estimated slope vs. True slope

9.768420 +/ 0.462170  17.006006

(That's 15.66 quantiles from true slope=17.006006)

True dP = 100.000, Eff. dP = 57.441+/2.718

NOW (time step 5 fs)
T P <V (A^3)> dV/<V>
44 797 36833.630 0.002426
44 897 36749.948 0.002431
compressibility from finite difference = 2.3046 atm^1
compressibility from fluctuation (at 797 atm) = 3.5693e05 atm^1
compressibility from fluctuation (at 897 atm) = 3.5772e05 atm^1
compressibility from fluctuation (average) = 3.5733e05 atm^1

Estimated slope vs. True slope

10.557018 +/ 0.301752  17.006006

(That's 21.37 quantiles from true slope=17.006006, FYI.)

True dP = 100.000, Eff. dP = 62.078+/1.774

~~~~~~~~~~~~~~~~
Michael Shirts
Associate Professor
michael.shirts@...780...
http://www.colorado.edu/lab/shirtsgroup/
Phone: (303) 7357860
Office: JSCBB C123
Department of Chemical and Biological Engineering
University of Colorado Boulder
On 7/19/17, 1:17 PM, "Michael R Shirts" <Michael.Shirts@...1657...> wrote:
Hi, all
I’m using MD to simulate FCC LennardJones. 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>/(((PhighPlow) (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.2565e05 atm^1
compressibility from fluctuation (at 797 atm) = 3.5237e05 atm^1
compressibility from fluctuation (at 897 atm) = 3.5108e05 atm^1
compressibility from fluctuation (average) = 3.5173e05 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.6792e05 bar^1
compressibility from fluctuation (at 797 bar) = 3.7367e05 bar^1
compressibility from fluctuation (at 897 bar) = 3.6073e05 bar^1
compressibility from fluctuation (average) = 3.6720e05 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:

Estimated slope vs. True slope

15.963630 +/ 0.345655  17.006006

(That's 3.02 quantiles from true slope=17.006006)

True dP = 100.000, Eff. dP = 93.871+/2.033

And for LAMMPS:

Estimated slope vs. True slope

9.768420 +/ 0.462170  17.006006

(That's 15.66 quantiles from true slope=17.006006)

True dP = 100.000, Eff. dP = 57.441+/2.718

Which seems problematic (GROMACS is a bit off, but much less so).
I’ve attached the lammps input files. To analyze them, I threw out the first 50000 steps for equilibration. I’ve attached the GROMACS mdp files (though not the coordinates and topologies because of size.
There are a few differences between the inputs, though they shouldn't affect the comparison I’m doing. I accidentally didn’t convert from bar to atm when going between GROMACS and LAMMPS runs, but did in the analysis. I don’t expect the average
volumes to be exactly the same between GROMACS and LAMMPS because of the way the cutoffs are handled, but they should be internally consistent within the program. The dimensions of the box are not quite the same between the two programs for various reasons,
but they are both roughly square, with 768 atoms. There are some difference in tau_t and tau_p, though they are all differences that shouldn’t really matter. Parameters should be the same – though again, it shouldn’t matter since we’re looking at internal
consistency. Here, GROMACS uses MartynaTuckermanTobiasKlein equations, and did not use mtk yes with lammps, though I would think the 1/V term should be negligible here (though I will test it).
Any insight would be useful!
Best,
~~~~~~~~~~~~~~~~
Michael Shirts
Associate Professor
michael.shirts@...780...
http://www.colorado.edu/lab/shirtsgroup/
Phone: (303) 7357860
Office: JSCBB C123
Department of Chemical and Biological Engineering
University of Colorado Boulder
