Fortran Simulation Lennard Jones
A large number of the examples simulate the Lennard-Jones liquid. Before discussing these in detail, we consider the different variants of the model that appear in the programs, and the state points used for testing.
For most of the examples, we use a cutoff of Rc=2.5σ. For MC programs the cut (but not shifted) potential (c) is used in the simulation, while for MD programs the cut-and-shifted potential (cs) is used. In all cases, an attempt is made to produce results without, and with, long-range corrections, the latter giving an estimate for the full potential (f). Differences between the different models are discussed in various places, see e.g.
Using their table V as a guide, we take the critical point to be roughly located at:
| LJ model | Tcrit | ρcrit | Pcrit |
|---|---|---|---|
| f: full, with LRC | 1.31 | 0.31 | 0.13 |
| c: cut (but not shifted) | 1.19 | 0.32 | 0.11 |
| cs: cut-and-shifted | 1.08 | 0.32 | 0.09 |
At any temperature below Tcrit, the liquid state is bounded below by the liquid-gas coexistence density, and using Tables II-IV of the same reference as a guide, we take the values for three example temperatures as
| LJ model | T | ρL | P |
|---|---|---|---|
| f: full, with LRC | 0.8 | 0.793 | 0.005 |
| c: cut (but not shifted) | 0.8 | 0.765 | 0.008 |
| cs: cut-and-shifted | 0.8 | 0.730 | 0.013 |
| LJ model | T | ρL | P |
|---|---|---|---|
| f: full, with LRC | 0.9 | 0.746 | 0.012 |
| c: cut (but not shifted) | 0.9 | 0.714 | 0.020 |
| cs: cut-and-shifted | 0.9 | 0.665 | 0.030 |
| LJ model | T | ρL | P |
|---|---|---|---|
| f: full, with LRC | 1.0 | 0.695 | 0.026 |
| c: cut (but not shifted) | 1.0 | 0.652 | 0.036 |
| cs: cut-and-shifted | 1.0 | 0.578 | 0.062 |
For Lennard-Jones, the state point (ρ,T)=(0.75,1.0) lies in the liquid region of the phase diagram for all these variants of the model.
Approximate equations of state for the full LJ potential have been presented by
- JK Johnson, JA Zollweg, KE Gubbins, Mol Phys, 78, 591 (1993),
- J Kolafa and I Nezbeda, Fluid Phase Equil, 100, 1 (1994),
- M Mecke, A Muller, J Winkelmann, J Vrabec, J Fischer, R Span, W Wagner, Int J Thermophys, 17, 391 (1996) and 19, 1493(E) (1998).
For testing our programs we used the more recent fitted equations of state presented by
- M Thol, G Rutkai, R Span, J Vrabec and R Lustig, Int J Thermophys, 36, 25 (2015),
- M Thol, G Rutkai, A Koester, R Lustig, R Span, J Vrabec, J Phys Chem Ref Data, 45, 023101 (2016),
for both the cut-and-shifted potential (denoted cs below), at Rc=2.5σ,
and the full potential (denoted f).
The formulae are implemented in the supplied program eos_lj.
For completeness, note that Thol et al also supply C++ programs, and tables of data,
in the Supplementary Information associated with their papers.
They are not responsible for our (Fortran) program!
First we look at the MD (and related) programs, which use the cut-and-shifted potential.
The first test of the MD codes is that energy, or the appropriate energy-like variable, is conserved. The following table uses runs of 10 blocks, each consisting of a number of steps to give 16 units of time per block with the indicated timestep (e.g. 1000×0.016, 2000×0.008 etc). We report the MSD values of the conserved variable for each program.
| δt | md_nve_lj |
md_nvt_lj |
md_npt_lj |
|---|---|---|---|
| 0.016 | 4.3284×10-6 | 4.0409×10-6 | 4.5875×10-6 |
| 0.008 | 1.8430×10-7 | 2.0036×10-7 | 2.4402×10-7 |
| 0.004 | 1.3121×10-8 | 1.4494×10-8 | 1.4956×10-8 |
| 0.002 | 1.2621×10-9 | 1.3526×10-9 | 1.5914×10-9 |
| 0.001 | 1.8530×10-10 | 2.1371×10-10 | 2.1005×10-10 |
Log-log plots show the expected dependence MSD ∝ δt4, hence RMSD ∝ δt2, except for some small deviations at the smallest timestep.
Now we compare EOS data
with typical test runs from our programs using default parameters, N=256, except where stated.
Note that E is the total internal energy per atom,
Cv is the constant-volume heat capacity per atom, and P is the pressure,
all including the ideal gas contributions.
The Smart Monte Carlo code smc_nvt_lj is included here since it uses the
cut-and-shifted potential which corresponds to the force calculation
(although it is not essential to do so).
Similarly, we include here the Brownian dynamics program bd_nvt_lj.
Numbers in parentheses (here and in the following tables) indicate errors in the last quoted digit, estimated from block averages. Results without error estimates are fixed (such as the temperature or density) or conserved.
| Source | ρ | T | E (cs) | P (cs) | Cv (cs) | E (f) | P (f) | Cv (f) |
|---|---|---|---|---|---|---|---|---|
| Thol et al (2015) (cs) | 0.75 | 1.00 | -2.9286 | 0.9897 | 2.2787 | |||
| Thol et al (2016) (f) | 0.75 | 1.00 | -3.7212 | 0.3996 | 2.2630 | |||
md_nvt_lj |
0.75 | 1.00 | -2.940(4) | 0.965(6) | 2.27(12) | -3.740(4) | 0.363(6) | 2.27(12) |
md_npt_lj§ |
0.7514(6) | 1.00 | -2.947(6) | 0.995(1) | -3.748(7) | 0.391(1) | ||
md_nve_lj |
0.75 | 1.0022(3) | -2.9289 | 0.987(2) | 2.24(1) | -3.7284 | 0.386(2) | |
md_nve_lj_omp |
0.75 | 1.0027(2) | -2.9278 | 0.986(2) | 2.28(1) | -3.7273 | 0.385(2) | |
md_nve_lj_vl |
0.75 | 1.0023(3) | -2.9278 | 0.992(2) | 2.24(1) | -3.7274 | 0.391(2) | |
md_nve_lj_ll‡ |
0.75 | 1.0010(1) | -2.9272 | 0.992(1) | 2.28(1) | -3.7268 | 0.391(1) | |
md_nvt_lj_ll‡ |
0.75 | 1.00 | -2.927(2) | 0.994(3) | 2.3(1) | -3.727(2) | 0.392(3) | 2.3(1) |
md_npt_lj_ll§‡ |
0.7501(4) | 1.00 | -2.931(4) | 0.991(1) | -3.731(4) | 0.389(1) | ||
smc_nvt_lj♯(a) |
0.75 | 1.00 | -2.9300(5) | 0.971(2) | 2.263(5) | -3.7296(5) | 0.369(2) | 2.270(5) |
smc_nvt_lj♯(b) |
0.75 | 1.00 | -2.928(2) | 0.99(1) | 2.26(2) | -3.728(2) | 0.39(1) | 2.27(2) |
smc_nvt_lj♯(c) |
0.75 | 1.00 | -2.930(3) | 0.98(2) | 2.26(3) | -3.729(3) | 0.38(2) | 2.27(3) |
bd_nvt_lj |
0.75 | 1.00 | -2.934(4) | 0.974(7) | 2.26(8) | -3.733(4) | 0.373(7) | 2.27(8) |
‡ Indicates a larger system size, N=864, needed to make the link-list method viable. Note that the speedup is not enormous for this system size, corresponding to 4×4×4 cells.
§ The constant-pressure simulations were run at P=0.99, the program default.
♯ The smc_nvt_lj program was tested (a) in default, single-particle-move, mode, with δt=0.1;
(b) in multi-particle mode, moving 100% of particles, with δt=0.02;
and (c) in multi-particle mode, moving 30% of particles, with δt=0.03.
These values give acceptance rates in the 45% – 55% range.
Results for md_lj_mts are not directly comparable,
because a larger cutoff (by default Rc=4.0σ) is used to illustrate the method.
The program was tested with N=400 (box length 8.1).
The usual state point is simulated: ρ=0.75 throughout.
No fitted EOS for the cs potential for this cutoff is available; obviously the estimates for the full potential
should be comparable with the values given above from Thol et al (2016).
Smallest timestep δt (called dt1 in the program)
and multiple-timestep-multipliers (the n_mts array) are given in column 2.
In all cases the run length was equal to 10 blocks of 20000 steps of length 0.005.
So, for example, 0.005 (142) means 10 blocks of 2500 steps of length 0.04,
each subdivided into 2 steps of length 0.02,
each subdivided again into 4 steps of length 0.005.
We also tried 0.002 (142) meaning 10 blocks of 6250 steps of length 0.016,
each subdivided into 2 steps of length 0.008,
each subdivided again into 4 steps of length 0.002.
The first line in the table below is from a run of md_nve_lj with the same system, for reference.
This example is just to illustrate the idea:
most of the test runs are actually slower, not faster, than md_nve_lj.
| Source | δt | T | E (cs) | P (cs) | Cv (cs) | E (f) | P (f) | E (MSD) |
|---|---|---|---|---|---|---|---|---|
md_nve_lj |
0.005 | 1.0038(1) | -3.5199 | 0.557(2) | 2.26(1) | -3.7157 | 0.410(2) | 1.7×10-8 |
md_lj_mts† |
0.005 (111) | 1.002(3) | -3.5199 | 0.58(1) | 2.4(1) | -3.7157 | 0.43(2) | 1.6×10-8 |
md_lj_mts‡ |
0.002 (142) | 1.0040(2) | -3.5196(2) | 0.559(1) | 2.26(1) | -3.7153(2) | 0.412(1) | 1.1×10-7 |
md_lj_mts§ |
0.005 (142) | 1.017(2) | -3.491(4) | 0.610(7) | 2.26(1) | -3.686(4) | 0.463(7) | 6.8×10-6 |
md_lj_mts¶ |
0.005 (142) | 1.0094(8) | -3.508(2) | 0.576(3) | 2.26(1) | -3.703(2) | 0.429(3) | 7.8×10-7 |
† All the timesteps the same length, as a check of the program book-keeping.
‡ Program default parameters; note the smaller timestep.
§ Program defaults, except for the timestep.
¶ Identical to § except that the switching length lambda is increased from 0.1 to 0.15.
Our MC programs use the cut (but not shifted) potential (which means that there is a delta correction, in the program, for the pressure). In this case, the value of Cv should be equal to the value for the full potential, since the energy LRC is independent of temperature. The Thol et al (2016) EOS for the full potential is used to predict results for the cut (but not shifted) potential (denoted c), at Rc=2.5σ, using the same LRC and delta corrections as in the MC codes. Once again, all values in the table include the ideal gas contribution. Except where indicated, tests are performed for N=256.
| Source | ρ | T | E (c) | P (c) | E (f) | P (f) | Cv (f) |
|---|---|---|---|---|---|---|---|
| Thol et al (2016) (f) | 0.75 | 1.00 | -3.3197 | 0.7008 | -3.7212 | 0.3996 | 2.2630 |
mc_nvt_lj |
0.75 | 1.00 | -3.332(1) | 0.651(3) | -3.734(1) | 0.350(3) | 2.28(1) |
mc_nvt_lj_re♯ |
0.75 | 1.00 | -3.332(1) | 0.648(2) | -3.734(1) | 0.347(2) | 2.258(4) |
mc_nvt_lj_ll‡ |
0.75 | 1.00 | -3.3320(3) | 0.659(2) | -3.7336(3) | 0.358(2) | 2.275(7) |
mc_npt_lj§ |
0.7501(2) | 1.00 | -3.331(1) | 0.666(2) | -3.733(1) | 0.364(2) | |
mc_npt_lj_llठ|
0.7499(2) | 1.00 | -3.330(1) | 0.663(1) | -3.732(1) | 0.362(1) | |
mc_zvt_lj¶ |
0.7504(4) | 1.00 | -3.333(3) | 0.668(4) | -3.735(3) | 0.366(4) | |
mc_zvt_lj_llঠ|
0.7506(3) | 1.00 | -3.335(2) | 0.671(3) | -3.736(2) | 0.370(3) |
‡ Indicates a larger system size, N=864 (or approximately so for mc_zvt_lj_ll).
Note that the linked lists do not give an enormous speedup for this system size,
which corresponds to 4×4×4 cells.
§ The constant pressure simulations were run at P=0.69, the program default.
The measured Cp (full) values were 5.28(7) for mc_npt_lj and 5.34(14) for mc_npt_lj_ll,
compared with Thol et al (2016) EOS giving 5.22.
The mc_npt_lj_ll program was run with non-default value db_max=0.015 to give a volume acceptance ratio around 20%.
¶ The grand canonical programs were run at activity z=0.0795, the program default value.
The Thol et al (2016) LRC-corrected value to give ρ=0.75 would be z=0.080627.
For mc_zvt_lj the box length was L=7σ; for mc_zvt_lj_ll L=10.5σ.
Acceptance rate of creation/destruction moves is quite small, at about 0.3%.
For other state points see below.
♯ The mc_nvt_lj_re program was run for four temperatures, see below for details.
The measured pressures P (c) are systematically a little low; this is particularly noticeable for the constant-pressure programs, where they might be expected to agree with the user-defined value of P. This reflects the approximate nature of the delta correction applied to the virial pressure, to account for the discontinuous potential at Rc. At the density ρ=0.75, with Rc=2.5, the pressure correction is Δ P≈-0.3, which is substantial. However, this estimate is based on the assumption that the pair distribution function g(Rc)=1. In fact, the choice Rc=2.5 is a poor one in this regard, lying near a local minimum where g(Rc)≈ 0.91 (an illustration of g(r) appears in the Fortran Analysis Structure page). Consequently the applied correction is slightly too large, and the resulting estimated pressure is systematically too low by ≈ 0.03. This serves as a reminder to always make clear what the cutoff is, and what corrections (for discontinuities or long-range interactions) have been applied.
In principle, there should be a delta correction for the configurational temperature. Long-range corrections to Tc are discussed by A Baranyai J Chem Phys, 112, 3964 (2000) and by A Lervik, O Wilhelmsen, TT Trinh, HR Nagel, J Chem Phys, 143, 114106 (2015).
Tests for the grand canonical MC program were initially conducted at a slightly lower density, very close to the liquid-vapour coexistence line (see Gibbs simulations below). A box length of 7σ was used, and creation/destruction acceptance ratios were around 1.5%. Comparison was made with the Thol et al (2016) equation of state, with corrections for the cutoff. The corresponding density is lower than the liquid coexistence density for the full potential, so there is no guarantee that the EOS will be accurate (it is only fitted in the single-phase regions of the full potential).
| Source | z | ρ | T | E (c) | P (c) | E (f) | P (f) |
|---|---|---|---|---|---|---|---|
| Thol et al (2016) (c) | 0.032 | 0.65325 | 1.0 | -2.7212 | 0.0457 | -3.0710 | -0.1828 |
mc_zvt_lj |
0.032 | 0.6532(5) | 1.0 | -2.728(3) | 0.0325(25) | -3.078(4) | -0.196(2) |
The program mc_gibbs_lj carries out Gibbs ensemble Monte Carlo,
and to test it we selected a temperature T=1.0,
which is below the critical point for the cut (but not shifted) LJ potential
(see tables above).
It was found convenient to start from a lower temperature,
with configurations at gas and liquid densities, with roughly equal numbers of particles,
and slowly work upwards in temperature, to equilibrate.
Note that the program expects two starting configurations: cnf1.inp and cnf2.inp.
The total number of atoms was fixed at NL+NG=512
and total volume VL+VG≈5514.
Exchanges of box identity are expected as the critical temperature is approached,
and so one should not place blind trust in the separate box averages reported by the program,
but refer to histograms of density, energy etc.,
illustrative examples of which appear here.
At T=1.0, however, these exchanges of box identity are expected to be infrequent, were not observed in the test runs, and the averages corresponded well to literature values for the coexistence parameters. The production run corresponded to default parameters in the program.
| Source | ρL | ρG | PL | PG | EL (c) | EG (c) |
|---|---|---|---|---|---|---|
| Trokhymchuk et al MC | 0.6542 | 0.0439 | 0.0336 | 0.0336 | ||
| Trokhymchuk et al MD | 0.6507 | 0.0500 | 0.0380 | 0.0380 | -2.713 ‡ | 1.047 ‡ |
mc_gibbs_lj |
0.653(1) | 0.050(1) | 0.031(2) | 0.038(1) | -2.731(5) | 1.049(9) |
‡ Indicates values for given ρ and T from the Thol et al (2016) EOS (f) with cutoff correction.
The small discrepancy between measured pressures in the two phases reflects the approximate nature of the delta correction for potential discontinuity, particularly in the liquid phase (see above). For a density ρ≈ 0.65 and Rc=2.5 the pressure correction is Δ P≈-0.23. However, this assumes g(Rc)=1, whereas actually g(Rc)≈ 0.95 at this density. Hence the correction is too large by approximately 0.01.
The mc_nvt_lj_re program uses MPI to handle communications between processes.
Here are some notes on the way the code is written.
We have only attempted to handle the most obvious errors at the start of the program, such as missing configuration files and incorrect user data, by closing down all the processes. A production code would take more care to handle exceptions during the run. Unhandled exceptions could possibly lead to processes hanging or becoming deadlocked, so you should be aware of the dangers in running this example.
In the program, all processes write to their standard output output_unit, but the default in MPI is
for all this output to be collated (in an undefined order) and written to a single channel. Testing
was carried out using Open MPI, which allows the program to be run with a command line which includes
an option for each process to write to separate files, similar to the following:
mpirun -np 4 -output-filename out ./mc_nvt_lj_re < mc.inp
whereby the output files are placed in subdirectories,
identified by process rank,
beneath the specified directory out
(in earlier versions of Open MPI, standard output would be directed to files named out##, the ## part being determined by the process rank).
If your implementation does not have this option, you should edit the code to explicitly open a file for
standard output, with a process-rank-dependent name, and associate the output_unit with it.
The mc_nvt_lj_re program conducts runs at several temperatures: four were used in testing.
The default program values include T=1.0, which is reported above, and here is the complete set,
with expected values from the Thol et al (2016) equation of state (f) corrected for cutoff.
As usual the program employed the cut (but not shifted) potential.
All runs are for density ρ=0.75, N=256, as usual.
At the lowest temperature, the full-potential system would lie in the coexistence region,
and the estimated pressure is negative.
| Source | T | E (c) | P (c) | E (f) | P (f) | Cv (f) |
|---|---|---|---|---|---|---|
| Thol et al (2016) (f) | 0.8772 | -3.6001 | 0.1942 | -4.0017 | -0.1070 | 2.3081 |
mc_nvt_lj_re |
0.8772 | -3.613(1) | 0.140(2) | -4.014(1) | -0.161(2) | 2.31(1) |
| Thol et al (2016) (f) | 1.0000 | -3.3197 | 0.7008 | -3.7212 | 0.3996 | 2.2630 |
mc_nvt_lj_re |
1.0000 | -3.332(1) | 0.648(2) | -3.734(1) | 0.347(2) | 2.258(4) |
| Thol et al (2016) (f) | 1.1400 | -3.0055 | 1.2571 | -3.4070 | 0.9559 | 2.2278 |
mc_nvt_lj_re |
1.1400 | -3.016(1) | 1.212(2) | -3.417(1) | 0.911(2) | 2.233(4) |
| Thol et al (2016) (f) | 1.2996 | -2.6523 | 1.8667 | -3.0539 | 1.5655 | 2.1989 |
mc_nvt_lj_re |
1.2996 | -2.662(1) | 1.820(3) | -3.063(1) | 1.519(3) | 2.214(5) |
The above (default) temperatures are chosen to give swap acceptance ratios all fairly close to 20% here
(of course, the set of temperatures, and all other run parameters, may be chosen by the user in a
namelist contained in the input file).
It should be noted that process m reports the swap acceptance ratio for exchanges with process m+1,
and the output file for the process with highest rank will report a zero swap ratio.
The programs md_nvt_lj_le and md_nvt_lj_llle are intended to illustrate:
the moving boundaries used in nonequilibrium shear flow simulations;
an algorithm for integrating the SLLOD equations of motion with constrained kinetic energy;
and an adapted link-cell method required to handle the modified boundaries.
These programs both use the short-ranged WCA Lennard-Jones potential,
in order to compare results with the following papers:
- G Pan, JF Ely, C McCabe, DJ Isbister, J Chem Phys, 122, 094114 (2005),
- KP Travis, DJ Searles, DJ Evans, Mol Phys, 95, 195 (1998).
Testing was performed at the state point used in those papers: ρ=0.8442, T=0.722. A system size N=256 was used (for such a short-range potential, this system size was also suitable for the link-cell program). The given program defaults, including a time step of 0.005, were used throughout, except for the strain rate which was varied.
| Strain rate | E | P | η | E | P | η |
|---|---|---|---|---|---|---|
| 0.04 | 1.8042(1) | 6.390(1) | 2.38(4) | 1.8039(3) | 6.389(2) | 2.31(4) |
| 0.16 | 1.8095(2) | 6.426(1) | 2.228(8) | 1.8098(2) | 6.427(1) | 2.23(1) |
| 0.64 | 1.8648(2) | 6.777(2) | 1.938(2) | 1.8646(2) | 6.776(1) | 1.935(2) |
In the table above, for each strain rate,
the results in columns 2-4 come from md_nvt_lj_le
and those in columns 5-7 from md_nvt_lj_llle
(essentially identical, but roughly twice as fast for N=256).
In all cases the kinetic energy was conserved very accurately by the algorithm.
The results, particularly the increase in E and P,
and the decrease in shear viscosity η,
as the strain rate increases,
are in good agreement with the above papers.
Incidentally, at the highest strain rate 0.64,
the configurational temperature is systematically about 1% lower
than the (constrained) kinetic temperature.