To investigate the structural, diffusive and thermodynamic properties of model systems in order to explore links between these properties. Particular emphasis is placed on uncovering the nature of diffusion kurtosis [1] in model systems.

MD simulations are traditionally performed under conditions of constant number of particles, volume and energy (the $$$N,V,E$$$-ensemble or micro canonical ensemble). The model we consider is that of a diffusing species (e.g. water) which is the solvent adsorbed in an immobile network, formed by spheres. The solvent particles interact with each other and the hard spheres. The number of solvent particles (species 1) is $$$N_{1}$$$, the number of solute particles (species 2) is $$$N_{2}$$$, the simulation box has side $$$L$$$ and volume $$$V=L^{3}$$$ and $$$L=\gamma R$$$, where $$$\gamma$$$ determines the size of the box. The colloid particles possess radius $$$R$$$ and the solvent particles are point-like. The number density of solvent particles is given by $$$\rho_{1}=\frac{N_{1}}{V}$$$ and the number of solute (colloid) particles is $$$\rho_{2}=\frac{N_{2}}{V}$$$. The reduced temperature is $$$T^{*}=1/\beta\epsilon$$$, where $$$\epsilon$$$ is the energy and $$$\beta=1/k_{B}T$$$ where $$$k_{B}$$$ is the Boltzmann constant and $$$T$$$ is the absolute temperature. We define $$$\sigma=\alpha R$$$ to be the length scale on which the solvent particles interact where $$$\alpha$$$ represents an effective measure of the size of the solvent to the size of the colloid. The solvent particles interact via the Lennard-Jones potential, given by $$u_{\textbf{LJ}}(r)= 4\epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left(\frac{\sigma}{r} \right)^{6} \right],$$ where $$$r$$$ is the distance between two interacting particles. It is a popular choice of potential in soft matter physics incorporating attractive ($$$(1/r)^{6}$$$ term) and repulsive ($$$(1/r)^{12}$$$ term) interactions. The resulting inter-particle force is $$F_{LJ}(r)=-\frac{du}{dr}=24\cdot\frac{\epsilon}{\sigma} \left[2 \left( \frac{\sigma}{r} \right)^{13} - \left(\frac{\sigma}{r} \right)^{7} \right].$$ The colloid particles are fixed in position, whereas the solvent particles are free to diffuse. This is realistic since over the timescale of the experiment ($$$T_{R}$$$, say) the colloid particles remain effectively fixed compared with the solvent particles. The inter-particle potential $$$v$$$ between colloid and solvent particles is infinite if the particles overlap and zero otherwise [2]. For numerical reasons we rather choose the potential:$${v(r)=\begin{cases}(\sigma/r)^{36}&\text{if particles overlap ($r \leq R$)} \\0&\text{otherwise ($r > R$), } \end{cases}}$$ where $$$r$$$ is the inter-particle distance and $$$R$$$ is the colloid radius. This represents a steep repulsion. The simulations (using C++) are performed by initiating the solvent particles at their natural crystallisation sites: a face-centred cubic lattice in this case. The colloid particles are then inserted randomly in the gaps. The solvent particles are assigned random velocities, subject to the conservation of total linear momentum, $$$\mathcal{P}$$$. The velocity Verlet algorithm [3] is used to propagate the system. Thermodynamic relationships are used to calculate quantities of interest. Periodic boundary conditions are used to approximate an infinite system. This means that the particle positions in all 27 "image" boxes must be considered, and interactions between a particle and the nearest image of another particle must be used for calculations.

Fig. 1 displays the instantaneous pressure and (negative) internal energy for a system of pure solvent particles. The conservation (aside from fluctuations) of these properties is an important check to test that the MD simulation is working. Since the particles are initiated with random velocities and on a lattice, which does not conform to the equilibrium state, the system takes some time to relax from an "equilibration" regime to a "production" regime. Information from the production regime is used to calculate properties of the system, including the pressure, energy, radial distribution function, diffusion and kurtosis.

In Fig. 2, a binary mixture of $$$N_{1} = 108$$$ solvent particles and $$$N_{2} = 49$$$ colloid particles is shown. The colloid particles remain spatially fixed during the simulation but the solvent particles are free to diffuse. The simulation has been terminated after 200 time steps so that the nature of the paths may be visualised before the trajectories fill too much space.

Fig. 3 displays the diffusion coefficient for a pure system with $$$N_{1} = 500$$$ and $$$T^{*} = 1.542$$$.

Fig. 4 displays the kurtosis for the same system as in Fig. 3. Note than in the long time limit, the kurtosis tends to zero, as expected for a pure system.

We have demonstrated that MD simulations may be used to investigate the diffusive, structural and thermodynamic properties of model systems. Some of these quantities (diffusion coefficient, kurtosis) are accessible by diffusion-weighted MRI. In addition to these preliminary investigations, future work involves performing more simulations to allow the examination of potential relationships between these properties, in particular the pressure and diffusion parameters.