-
Notifications
You must be signed in to change notification settings - Fork 105
Python Simulation Lennard Jones
State points for the Lennard-Jones potential, in its cut (but not shifted) form, denoted (c), cut-and-shifted form, denoted (cs), and full form (f), are discussed in the Fortran example. Except where otherwise stated, we use our default liquid state point (ρ,T)=(0.75,1.0) for testing, with N=256 atoms, Rc=2.5σ, and compare with the same equations of state due to Thol et al
- 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).
We provide a Python version of the equation-of-state program eos_lj.py
based on
the formulae in those papers.
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 (Python) program!
The Python codes run much more slowly than the Fortran ones,
and so typically our default parameters carry out shorter runs
(e.g. 10 blocks of 1000 steps rather than 10 blocks of 20000 steps).
For both MC and MD examples,
we supply two versions of the code in the main module:
a slow version built around the standard Python loops,
which may be directly compared with the Fortran examples;
and a fast version which uses NumPy library routines to replace the inner loop.
This is tantamount to vectorizing the calculation in a very simple way.
(We do not claim to have chosen the fastest NumPy method in our Python simulation codes;
for examples of other approaches to the same problem, see the pair_distribution.py
program).
The variable fast
controls which version is used, and by default is set to True
,
at the start of the main module; the user may override this by editing this statement,
but the program will then run (yet another) order of magnitude more slowly,
so it may be necessary to reduce the run length still further.
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.1613×10-6 | 4.3865×10-6 | 5.4746×10-6 |
0.008 | 1.8896×10-7 | 2.0417×10-7 | 2.5920×10-7 |
0.004 | 1.3705×10-8 | 1.3769×10-8 | 2.0192×10-8 |
0.002 | 1.3866×10-9 | 1.1486×10-9 | 1.5887×10-9 |
0.001 | 1.2650×10-10 | 1.8297×10-10 | 2.3865×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 of our programs. The results in the following table use the same run lengths as for the Fortran examples (i.e. longer than the default value in the program), and default parameters otherwise. Results for Smart Monte Carlo and Brownian Dynamics are included, as they use the same cut-and-shifted potential.
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_nve_lj.py |
0.75 | 1.0027(4) | -2.9278 | 0.988(3) | 2.24(1) | -3.7274 | 0.387(3) | |
md_nvt_lj.py |
0.75 | 1.00 | -2.937(3) | 0.975(4) | 2.1(1) | -3.737(3) | 0.374(4) | 2.1(1) |
md_npt_lj.py § |
0.7509(5) | 1.00 | -2.942(5) | 0.994(1) | -3.743(6) | 0.3908(6) | ||
md_nve_lj.py ‡ |
0.75 | 1.0006(1) | -2.9283(1) | 0.9903(7) | 2.26(1) | -3.7279(1) | 0.3887(7) | |
md_nvt_lj.py ‡ |
0.75 | 1.00 | -2.934(2) | 0.982(3) | 2.4(1) | -3.734(2) | 0.381(3) | 2.4(1) |
md_npt_lj.py §‡ |
0.7505(4) | 1.00 | -2.936(4) | 0.9920(7) | -3.736(5) | 0.3896(2) | ||
smc_nvt_lj ♯(a) |
0.75 | 1.00 | -2.929(1) | 0.979(3) | 2.256(4) | -3.729(1) | 0.378(3) | 2.264(4) |
smc_nvt_lj ♯(b) |
0.75 | 1.00 | -2.932(2) | 0.95(1) | 2.28(3) | -3.732(2) | 0.35(1) | 2.29(3) |
smc_nvt_lj ♯(c) |
0.75 | 1.00 | -2.934(3) | 0.94(1) | 2.24(2) | -3.733(3) | 0.34(1) | 2.24(2) |
bd_nvt_lj.py |
0.75 | 1.00 | -2.931(3) | 0.976(5) | 2.2(1) | -3.731(3) | 0.374(5) | 2.2(1) |
§ The constant-pressure simulation was 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.
‡ Indicates a larger system size, N=864, and a neighbour-list method:
the main program is edited to import routines from md_lj_ll_module.py
rather than md_lj_module.py
.
A cell structure is used, similar to the Fortran example.
Linked lists are not used (despite the module name!), as there is no advantage to this in Python.
The resulting code is not especially efficient for this system size, with 4×4×4 cells,
and in a practical application when performance is an issue
Python/NumPy would probably not be used for the force routine.
However, the module illustrates the basic idea of this approach.
The results in the following table use the same run lengths as for the Fortran examples, and default parameters otherwise.
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.py |
0.75 | 1.00 | -3.3315(7) | 0.653(4) | -3.7331(7) | 0.352(4) | 2.274(9) |
mc_nvt_lj_re ♯ |
0.75 | 1.00 | -3.3314(4) | 0.652(2) | -3.7330(4) | 0.351(2) | 2.272(5) |
mc_npt_lj.py |
0.7510(6) | 1.00 | -3.339(4) | 0.69 | -3.741(5) | 0.360(2) | |
mc_zvt_lj.py |
0.7497(3) | 1.00 | -3.329(2) | 0.656(4) | -3.731(2) | 0.355(4) |
♯ 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 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 Python 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.
The program mc_gibbs_lj.py
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 ready-equilibrated configurations
generated by the corresponding Fortran code,
and these in turn had been produced by starting from lower temperatures,
and slowly raising T.
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, except that the run length was taken to be 10 blocks of 10000 steps, as for the Fortran example.
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.py |
0.653(1) | 0.050(1) | 0.028(1) | 0.038(1) | -2.729(5) | 1.053(7) |
‡ 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.py
program uses MPI to handle communications between processes.
Here are some notes on the way the code is written.
We have used the mpi4py package of Python bindings to MPI. To run this example you will have to install this package, with a compatible implementation of MPI. We have tested it with Open MPI. Other Python packages which interface with MPI and other implementations of MPI itself are, of course, available.
We have adopted the simple approach of communicating single Python objects
using the supplied methods with lowercase names such as bcast
, send
, recv
, allreduce
,
which employ pickle behind the scenes.
NumPy arrays are exchanged as buffers, using the uppercase-named method Sendrecv_replace
,
with automatic detection of data type and length.
This should be faster, but we have not made much effort to optimize the speed.
The MPI Init
and Finalize
functions are not explicitly called:
mpi4py invokes them automatically when the MPI module is loaded, and when the Python process exits, respectively.
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, 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:
mpiexec -n 4 -output-filename out ./mc_nvt_lj_re.py < 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 channel with it.
The mc_nvt_lj_re.py
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;
default parameters are used, except for the run lengths,
which are set to 10 blocks of 10000 steps, as per the Fortran version.
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.6130(6) | 0.141(3) | -4.0146(6) | -0.160(3) | 2.310(7) |
Thol et al (2016) (f) | 1.0000 | -3.3197 | 0.7008 | -3.7212 | 0.3996 | 2.2630 |
mc_nvt_lj_re |
1.0000 | -3.3314(4) | 0.652(2) | -3.7330(4) | 0.351(2) | 2.272(5) |
Thol et al (2016) (f) | 1.1400 | -3.0055 | 1.2571 | -3.4070 | 0.9559 | 2.2278 |
mc_nvt_lj_re |
1.1400 | -3.0155(4) | 1.211(2) | -3.4171(4) | 0.910(2) | 2.224(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.6620(7) | 1.819(3) | -3.0636(7) | 1.518(3) | 2.206(6) |
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
JSON list 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 program md_nvt_lj_le.py
is intended to illustrate
the moving boundaries used in nonequilibrium shear flow simulations and
an algorithm for integrating the SLLOD equations of motion with constrained kinetic energy.
The program uses 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.
The given program defaults, including a time step of 0.005,
were used to give the results in the table below,
except that the run length was increased to 10×100000 steps,
the same as for the Fortran example,
and the strain rate was varied as shown.
In the table below, for each strain rate,
the results in columns 2-4 come from md_nvt_lj_le.py
.
Strain rate | E | P | η | E | P | η |
---|---|---|---|---|---|---|
0.04 | 1.8040(2) | 6.389(1) | 2.31(3) | 1.8035(5) | 6.387(3) | 2.4(1) |
0.16 | 1.8099(3) | 6.428(2) | 2.227(9) | 1.8087(4) | 6.421(3) | 2.19(2) |
0.64 | 1.8648(2) | 6.777(1) | 1.940(2) | 1.8645(4) | 6.775(3) | 1.933(6) |
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, and with the results of the Fortran programs. Incidentally, at the highest strain rate 0.64, the configurational temperature is systematically about 1% lower than the (constrained) kinetic temperature.
The results in columns 5-7 come from the same program,
edited to import routines from md_lj_llle_module.py
rather than md_lj_le_module.py
.
This version uses neighbour lists, in the same way as md_lj_ll_module.py
.
The system size is unchanged;
however, for this system the program is very inefficient (compared with the simpler version without lists)
and consequently the runs used for testing are only 10×20000 steps.
The inefficiency probably results from the very small cells and hence short NumPy arrays,
which incur significant overheads.
It is worth repeating that in a practical application,
when performance is an issue,
Python/NumPy would probably not be used for the force routine.