Python Simulation Hard Particles - Allen-Tildesley/examples GitHub Wiki

Hard-particle programs

The programs mc_nvt_hs.py and md_nve_hs.py illustrate, respectively, the simplest MC and MD methods for the basic hard-sphere model. The temperature is not important in the first case: a factor kT is used to normalize the energies. The energy, in the second case, is identical with the (exactly conserved) kinetic energy, and hence closely related to the temperature. Equations of state for this model have been reported many times. Here we refer to some fairly recent, useful, sources of data and/or fitted equations

The paper of Kolafa et al (2004) is particularly careful to discuss corrections due to different ensembles and system size. Here we just present the raw results for a small system, N=256; programs are run with default parameters, except that the test runs were longer than default: 10 blocks of 10000 steps, same as for the Fortran examples. Starting fcc lattice configurations may be prepared using initialize.py in the usual way. The EOS is taken from the Hansen-Goos (2016) paper, and a program to evaluate it may be found in eos_hs.py.

ρ P (EOS) P mc_nvt_hs.py P md_nve_hs.py ρ mc_npt_hs.py
0.50 1.6347 1.626(2) 1.633(1) 0.499(2)
0.55 2.0574 2.053(3) 2.056(1) 0.552(2)
0.60 2.5769 2.571(3) 2.573(1) 0.600(2)
0.65 3.2171 3.207(6) 3.215(1) 0.642(1)
0.70 4.0087 3.994(6) 4.005(2) 0.702(2)
0.75 4.9910 4.959(6) 4.984(1) 0.752(3)

We must remember that P is calculated by a box-scaling method in the NVT simulation, which may introduce a small systematic error. This can be reduced by reducing the scaling factor, at the expense of worsening the statistics. We also provide a program mc_npt_hs.py to illustrate the constant-NPT method. For the averages of ρ reported above, the input pressure was that given by the corresponding EOS entry. With default parameters, volume move acceptance ratio was nearly 5% at the highest pressure, and around 11% at the lowest pressure studied here.

We also provide two programs to simulate the hard spherocylinder model, of cylinder length L and diameter D: mc_npt_sc.py and mc_nvt_sc.py. Configurations may be prepared as described in the Fortran example. Test runs were performed using 10 blocks of 10000 steps (as for the Fortran examples); the program default is 10×1000. For L=5, N=256 (a very small system, not recommended for serious work) we compare with the results of

See the Fortran example for comments about units, and other literature values.

P vmol ρ vmol P ρ S ρ S P S
(M) (M) (M) (M) (M) mc_npt_sc mc_npt_sc mc_nvt_sc mc_nvt_sc
2.53 0.310 0.568 0.070 0.041 0.0693(3) 0.096(8) 0.577(2) 0.089(5)
3.63 0.352 0.816 0.079 0.053 0.0794(2) 0.119(8) 0.814(2) 0.081(8)
4.89 0.397 1.099 0.089 0.136 0.0898(2) 0.32(1) 1.098(3) 0.24(2)
5.05 0.400 1.135 0.090 0.170 0.0901(2) 0.191(8) 1.141(5) 0.15(2)‡
5.40 0.419 1.213 0.094 0.574 0.0948(2) 0.57(1) 1.224(4) 0.54(1)
5.80 0.436 1.303 0.098 0.714 0.0991(2) 0.769(8) 1.281(6) 0.791(3)
6.20 0.448 1.393 0.101 0.754 0.1024(2) 0.834(6) 1.381(5) 0.789(5)

The mc_npt_sc runs use pressures from column 3 above; the mc_nvt_sc runs are at densities taken from column 4. At the highest pressure, using default parameters, move acceptance ratio was around 30%, and volume acceptance ratio around 10%. These values rose to 50% and 15% respectively at the lowest pressure. The ‡ reminds us that results for these run lengths, near the isotropic-nematic transition, where very slow evolution of the nematic order parameter would be observed, are unreliable. Also the system size is about 25% that used by McGrother, which has a direct effect on the measured nematic order parameter. With these caveats in mind, agreement between the two programs, and with the results of McGrother et al (1996), is reasonable.

⚠️ **GitHub.com Fallback** ⚠️