Fortran Simulation Hard Particles

Hard-particle programs

The programs mc_nvt_hs and md_nve_hs 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

• H Hansen-Goos, J Chem Phys, 144, 164506 (2016),
• MN Bannerman, L Lue, LV Woodcock, J Chem Phys, 132, 084507 (2010),
• J Kolafa, S Labik, A Malijevsky, Phys Chem Chem Phys, 6, 2335 (2004).

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. Starting fcc lattice configurations may be prepared using initialize 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.f90.

ρ	P (EOS)	P mc_nvt_hs	P md_nve_hs	ρ mc_npt_hs
0.50	1.6347	1.634(2)	1.632(1)	0.502(2)
0.55	2.0574	2.051(3)	2.055(1)	0.553(3)
0.60	2.5769	2.573(4)	2.573(1)	0.600(3)
0.65	3.2171	3.210(8)	3.215(2)	0.651(3)
0.70	4.0087	3.996(8)	4.005(3)	0.700(2)
0.75	4.9910	4.960(7)	4.985(4)	0.749(2)

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 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 and mc_nvt_sc. In preparing configurations for these programs, one must not allow overlap between the hard particles. A simple approach is to run initialize with molecules="linear", lattice=.false., and to request a very low density. For the default length=5 spherocylinders (L=5, D=1) a value of density=0.05 is suitable. Then, the constant-pressure program may be used to compress the system to higher density. This is a fairly slow process, requiring the density ρ and nematic order parameter S to be carefully monitored. Once suitable high-density state points have been prepared, a configuration at a precisely specified density, for use in the constant-volume program, may be obtained by a small expansion (using the adjust program).

For testing we use N=256; such a small system is not recommended for serious work, but allows us to explore up to ρ=0.148 (box length 12 D, twice the interaction range of 6 D) which is sufficient for our purposes. Equations of state from MC simulations are presented in two papers

• SC McGrother, DC Williamson, G Jackson, J Chem Phys, 104, 6755 (1996),
• PG Bolhuis, D Frenkel, J Chem Phys, 106, 666 (1997).

In making comparisons, care must be taken with the units. McGrother et al (1996) quote densities in the form of the packing fraction η=ρ v_mol and pressures as P v_mol, where v_mol is the molecular volume. We translate selected values from their Table V (denoted (M) below) into our reduced units based on D=1 below; for L/D=5, v_mol=4.4506. (Bolhuis and Frenkel (1997) define reduced densities relative to the density of closest packing of spherocylinders, while reporting pressures the same way as McGrother et al. We do not use the Bolhuis-Frenkel results below.)

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.0698(2)	0.081(7)	0.579(2)	0.073(6)
3.63	0.352	0.816	0.079	0.053	0.0799(2)	0.098(6)	0.810(3)	0.098(6)
4.89	0.397	1.099	0.089	0.136	0.0895(2)	0.25(1)	1.126(5)	0.23(1)
5.05	0.400	1.135	0.090	0.170	0.0903(2)	0.32(3)	1.136(3)	0.149(7)‡
5.05	0.400	1.135	0.090	0.170	0.0923(3)	0.502(7)	1.061(3)	0.592(8)‡
5.40	0.419	1.213	0.094	0.574	0.0966(3)	0.764(6)	1.193(6)	0.596(3)
5.80	0.436	1.303	0.098	0.714	0.0984(2)	0.783(8)	1.325(6)	0.751(6)
6.20	0.448	1.393	0.101	0.754	0.1021(3)	0.839(3)	1.406(6)	0.795(4)

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. Several successive runs (10 blocks of 10000 steps each) were undertaken for state points near the isotropic-nematic transition, where very slow evolution of the nematic order parameter could be observed. Considerable sluggishness is expected around this transition, giving irreproducible results, an example indicated by ‡ in the table. Much longer runs are needed to achieve consistency! 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, which apply mainly to the middle three state points in the table, agreement between the two programs, and with the results of McGrother et al (1996), is reasonable.