System properties

The energy in an MD simulation consists of the kinetic energy:

\[E_{kin} = \frac{1}{2} \sum_i m_i v_i^2\]

and the potential (configurational) energy:

\[E_{pot} = U_c = \sum_i \sum_{j<i} V \left(r_{ij}\right).\]

with both contributions interchanging during the course of the simulation. Both properties can be tracked and ensemble-averaged as static system properties, along with:

• Temperature

\[T = \frac{2 \langle E_{kin} \rangle}{3 k_B}\]

• Pressure

\[P = N k_B T - \frac{1}{3} \Big \langle \sum_i^{N} \vec{r}_i \cdot \vec{f}_i \Big \rangle\]

• Specific heat

\[\langle \delta \left(U_c \right) \rangle = \frac{3}{2} N k_B T \left(1 - \frac{3 N k_B}{2 C_v} \right)\]

noting that these values typically fluctuate during simulations and should be reported with error estimates.

The structure of a material can be described by a pair correlation function, the radial distribution function:

\[g \left( r \right) = \frac{\langle n \left( r \right) \rangle}{4 \pi \rho r^2 \Delta r} = \frac{V}{N^2} \Bigg \langle \sum_{i=1}^{N-1} \sum_{j \neq i}^{N} \delta \left(r - r_{ij}\right) \Bigg \rangle\]

which can be calculated by finding distances between atom pairs and keeping count of the numbers, \(n \left( r \right)\), in ‘bins’ of size \(\Delta r\). The radial distribution function can be related to the structure factor:

\[S \left( k \right) = 1 + 4 \pi \rho \int_0^{\infty} \frac{\sin \left(k r\right)}{k r} \left(g \left( r \right) - 1 \right) r^2 dr\]

which can also be obtained experimentally from X-ray diffraction.

MD simulations can also be used to find various dynamic properties, including mean squared displacements (MSDs):

\[\langle | \vec{r}_i \left(t\right) - \vec{r}_i \left(0\right)|^2 \rangle\]

and velocity autocorrelation functions (VAFs):

\[\langle \vec{v}_i \left(t\right) \cdot \vec{v}_i \left(0\right) \rangle\]

both of which can be related to the diffusivity of the material:

\[D = \frac{\langle | \vec{r}_i \left(t\right) - \vec{r}_i \left(0\right)|^2 \rangle}{6t} = \frac{1}{3} \int_0^{t} \langle \vec{v}_i \left(t\right) \cdot \vec{v}_i \left(0\right) \rangle dt\]

The shear viscosity can be found by using shear stress autocorrelation functions (SSAFs):

\[\mu_0 = \frac{V}{k_B T} \int_0^t \langle \sigma_{\alpha \beta} \left(t\right) \cdot \sigma_{\alpha \beta} \left(0\right) \rangle dt\]

based on off-diagonal components of the stress tensor \(\sigma_{\alpha \beta}\) (i.e. \(\alpha \neq \beta\)). This value is obtained at the limit of zero shear rate.

Alternatively, the relationship between shear stress and velocity gradient (shear rate):

\[\tau = \mu \frac{\partial u_x}{\partial y}\]

can be determined using non-equilibrium MD simulations to find the shear viscosity and the system rheology. These involve applying constant velocity gradients and measuring the resulting shear stresses, given as the relevant component of the pressure tensor (stress tensor), \(\tau = |\langle\sigma_{yx}\rangle|\). Simple Newtonian fluids produce a constant viscosity regardless of shear rate, while longer molecules and solutions often do not.

^ GO TO TOP ^