Dear Fellow Lammpers,
I have been computing the chemical potential at fixed volume and temperature of a simple Lennard Jones system (density around 0.75 and temperature 2 in Lennard Jones units) and comparing my estimates of the chemical potential depending on whether I use Particle mesh ewald for dispersion, or impose a simple cut-off at a range of cut-off values (RCUT = 3,4.5,6,7.5,9,11,13 and 14.5), and with published data where the estimates were obtained using Equation of state methods.
My real interest was to use PME – but I needed to be sure I was getting convergence between either approach – which should happen if RCUT is sufficiently long. Unfortunately, the PME estimates of the chemical potential are seriously wrong, while the simple approach converges to published values for RCUT >=6 (RCUT = 4.5 is already pretty close). Either I am making an error in my setup for PMME or there is a serious error in the code. Anyone come across something wierd – or am I making a common error? Some more details on the run conditions and the estimates of the chemical potential are below.
I made the obvious checks such as comparing the potential energy for the same coordinates using PME and LJ/CUT (with a huge rcut = 14) – they agree for the first signifcant 3 digits – hence I am pretty baffled why there is such a large difference for the corresponding chemical potentials. I should mention that to compute the chemical potential, I used effectively a form of thermodynamic perturbation where the sigma of the Lennard jones of inserted particles is increased from 0 to 1 in ten steps – maybe something in the grid representation of dispersion is tripping things up. At the bottom of this message are the results for a short NVE run to perform a sanity check - suggesting that all should be fine - but alas is not.
Thanks for any wisdom/suggestions,
Donal
I am using lammps 30Jul16 to estimate the chemical potential using machines at our national hpc centre…
neighbor 0.2 bin
neigh_modify delay 0 every 5 check yes one 12000 page 200000
variable rrcut index 4.5
pair_style lj/long/coul/long long off {rrcut}
pair_modify mix geometric
pair_coeff 1 1 1.0 1.0 {rrcut}
pair_coeff 2 2 1.0 1.0 ${rrcut}
kspace_style pppm/disp 0.001
kspace_modify mesh 10 10 10 order 5 gewald/disp 0.65 mix/disp geom
To compare without using particle mesh ewald, I simple set
pair_style lj/long/coul/long off off ${rrcut}
and increase rrcut to see when I get convergence in the estimates for the chemical potential.
To compute the chemical potential, every one hundred time steps, I compute the potential energy for ten intermediate values of a perturbation parameter (the sigma of an inserted particle) (using run 0 to force the evaluation for each instance)
Chemical Potential Estimates
density 0.745283138 temp 2 (3200 atoms)
EOS 3.080689 (published result – Johnson et al)
MBAR (used to estimate chemical potential)
PME 4.1592415336058899 error bar 0.062103390104848596
RCUT3.0 3.349640795197935 error bar0.061616726883664964
RCUT4.5 3.1107173736649827 error bar 0.062154241611203975
RCUT6.0 3.0858444532817799 error bar 0.062259115856097377
RCUT7.5 2.9893428943105214] error bar 0.062052339461932927
RCUT8.5 3.0160919804911472 error bar 0.062164063898934224
RCUT10.0 3.0045759611025389 error bar 0.061879346043870874
RCUT11.5 3.0216386890807065 error bar 0.06196124877978211
RCUT13.0 3.0103565714031006 error bar 0.0621380232603004
RCUT14.5 3.064600538159425 error bar 0.062421311951306874
Comarison of output for NVE (a few steps) using same initial condition (PME versus LJ/CUT with rcut = 14) run on my laptop as a quick check.
PME RCUT = 4.5
Dispersion G vector (1/distance)= 0.65
Dispersion grid = 18 18 18
Dispersion stencil order = 5
LAMMPS (17 Dec 2016)
Reading restart file …
restart file = 30 Jul 2016, LAMMPS = 17 Dec 2016
WARNING: Restart file used different # of processors (…/read_restart.cpp:717)
orthogonal box = (8.51011 5.18889 -0.897307) to (24.7635 21.4422 15.356)
2 by 2 by 1 MPI processor grid
3200 atoms
Finding 1-2 1-3 1-4 neighbors …
Special bond factors lj: 0 0 0
Special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
3200 atoms in group real
Resetting global fix info from restart file:
fix style: nvt, fix ID: 1
PPPMDisp initialization …
WARNING: Charges are set, but coulombic solver is not used (…/pppm_disp.cpp:328)
Dispersion G vector (1/distance)= 0.65
Dispersion grid = 18 18 18
Dispersion stencil order = 5
Dispersion estimated absolute RMS force accuracy = 0.000606772
Dispersion estimated absolute real space RMS force accuracy = 9.29222e-05
Dispersion estimated absolute kspace RMS force accuracy = 0.000599615
Dispersion estimated relative force accuracy = 0.000606772
using double precision FFTs
3d grid and FFT values/proc dispersion = 4508 1620
All restart file global fix info was re-assigned
Neighbor list info …
update every 5 steps, delay 0 steps, check yes
max neighbors/atom: 12000, page size: 200000
master list distance cutoff = 4.7
ghost atom cutoff = 4.7
binsize = 2.35, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/long/coul/long, perpetual
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Setting up Verlet run …
Unit style : lj
Current step : 300000
Time step : 0.00471524
Memory usage per processor = 6.91285 Mbytes
Step Press Temp TotEng PotEng Volume KinEng
300000 3.9940242 1.9943461 -1.5160062 -4.5065904 4293.6703 2.9905843
300001 4.0111161 1.9913578 -1.5156816 -4.5017848 4293.6703 2.9861033
300002 4.0157406 1.9910637 -1.5154572 -4.5011194 4293.6703 2.9856622
300003 4.008061 1.9934185 -1.5153514 -4.5045447 4293.6703 2.9891933
300004 3.9932769 1.9972575 -1.5152811 -4.5102311 4293.6703 2.99495
300005 3.9784125 2.0009975 -1.5151582 -4.5157166 4293.6703 3.0005584
300006 3.9687031 2.0034437 -1.5149635 -4.5191898 4293.6703 3.0042264
300007 3.9664534 2.0040618 -1.514718 -4.5198713 4293.6703 3.0051533
300008 3.9715514 2.0028757 -1.5144294 -4.5178042 4293.6703 3.0033747
300009 3.981785 2.0003872 -1.5140994 -4.5137425 4293.6703 2.9996431
300010 3.9930874 1.997514 -1.513751 -4.5090856 4293.6703 2.9953346