Rep:Mod:sfs114

Theory

Numerical Integration

The classical solution for the position at time $t$ compares well with the velocity-Verlet solution:

A plot of x(t) against t comparing analytical and velocity-Verlet solutions

The total energy for the oscillator varies as shown:

A plot of energy against time

An approximate linear fit has been performed on the maxima of the error of the calculations; the absolute difference between classical and velocity-Verlet solutions. Iterations of using previous results causes error to propagate and increase.

A plot of absolute error, with a linear function fitted to the error maxima

The smaller the timestep, the smaller fluctuations in total energy. Calculations over larger timesteps causes a greater error, as particles could end up too close together and face extremely large forces, for example. It is important to monitor the total energy of a physical system to ensure energy is conserved, however infinitesimal timesteps greatly increase time needed to run simulations. Larger timesteps allow a longer length of time to be simulated. A timestep of 0.028s allows energy fluctuations to be as low as $\pm 1\mathrm{\%}$ of the average value.

Energy oscillation with a timestep of 0.028s

Energy oscillation with a timestep of 0.5s

Atomic Forces

For a single Lennard-Jones interaction, $\varphi (r)=4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})$.

• When potential energy is 0, φ(r) = 0 and r=r₀

If $\varphi =0,(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})=0$,

$\frac{{\sigma }^{12}}{r^{12}}=\frac{{\sigma }^6}{r^6}$

$r_0=\sigma$

• The force is given by $F=\frac{d\varphi }{dr}=4ϵ(-12{\sigma }^{12}{r}^{-13}+6{\sigma }^6{r}^{-7})$ and at a potential energy of 0, $F=\frac{-24ϵ}{\sigma }$.

• At equilibrium separation $r_{eq}$, $\frac{d\varphi }{dr}=0$.

$0=4ϵ(\frac{-12{\sigma }^{12}}{r^{13}}+\frac{6{\sigma }^6}{r^7})$

$(\frac{12{\sigma }^{12}}{r^{13}}=\frac{6{\sigma }^6}{r^7})$

$r_{eq}=2^{\frac{1}{6}}\sigma$

• At equilibrium separation, the well depth:

$\varphi (r_{eq})=\varphi (2^{\frac{1}{6}}\sigma )=-ϵ$

• $\int \varphi (r)dr=4ϵ[\frac{-{\sigma }^{12}}{11r^{11}}+\frac{{\sigma }^6}{5r^5}]+c$ and given that $\sigma =ϵ=1.0$ so $\int \varphi (r)dr=4[\frac{-1}{11r^{11}}+\frac{1}{5r^5}]+c$

• $\underset{2\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi (r)dr=4{[\frac{-1}{11r^{11}}+\frac{1}{5r^5}]}_2^{\mathrm{\infty }}=-4(\frac{1}{5(2{)}^5}-\frac{1}{11(2{)}^{11}})=-0.02482$
• $\underset{2.5\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi (r)dr=4{[\frac{-1}{11r^{11}}+\frac{1}{5r^5}]}_{2.5}^{\mathrm{\infty }}=-4(\frac{1}{5(2.5{)}^5}-\frac{1}{11(2.5{)}^{11}})=-0.008177$
• $\underset{3\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi (r)dr=4{[\frac{-1}{11r^{11}}+\frac{1}{5r^5}]}_3^{\mathrm{\infty }}=-4(\frac{1}{5(3{)}^5}-\frac{1}{11(3{)}^{11}})=-0.00329$

Periodic Boundary Conditions

As $pV=N{k}_{B}T$, the number of water molecules in 1 mL of water is approximately $2.46x10^{19}$ and 10000 molecules takes up an approximate volume of $4.06x10^{-22}{m}^3$.

In a simulation box which runs from (0,0,0) to (1,1,1), an atom that starts at (0.5,0.5,0.5) and moves along vector (0.7,0.6,0.2), will end up at (0.2,0.1,0.7) once periodic boundary conditions have been applied.

Reduced Units

The LJ parameters for Argon are: $\sigma =0.34nm,\frac{ϵ}{k_B}=120K$.

• $r=\sigma \cdot r*=0.34\cdot 3.2=1.088nm$.
• $WellDepth=\varphi (r_{eq})=-ϵ=-k_B\cdot 120=-1.656^{-21}J$

$-1.656^{-21}\cdot N_A=-997.4=-0.997kJmo{l}^{-1}$.

• $T=\frac{ϵ\cdot T*}{k_B}=120\cdot 1.5=180K$.

Equilibriation

Creating the Simulation Box

If two atoms are generated too close together, the LJ potential shows that the potential between the two would be infinitely large, making force simulations between these two atoms too large to realistically simulate. The LJ cutoff also ensures that LJ potentials are only calculated for atoms that are near enough, and not every other atom in the infinitely repeating lattice, which would greatly increase simulation run time.

A lattice spacing of 1.07722 corresponds to a lattice number density of $\frac{1}{1.07722^3}=0.8$ for a simple cubic lattice. A face centred cubic lattice has 4 lattice points per cell, and thus would require a lattice spacing of $(\frac{4}{1.2}{)}^{\frac{1}{3}}=1.4938$. A 10x10x10 box would contain 1000 unit cells, and 4000 lattice points, so the create_atoms command for such a lattice would create 4000 atoms.

Setting the Properties of the Atoms

The command mass 1 1.0 assigns all atoms of type 1 a mass of 1.0. The command pair_style lj/cut 3.0 defines the cutoff distance between atoms that have a potential between them to be 3.0 (ie. the simulation does not run for atoms farther apart or closer than this distance). The command pair_coeff ** 1.0 1.0 specifically defines the pairwise force field coefficients for multiple pairs atoms.

The velocity-Verlet algorithm is the numerical integration method that will be used if $x_i(0)$ and $v_i(0)$ are defined.

Running the Simulation

Calling upon variables, instead of assigning numbers, makes it much easier to change these variables for every simulation that is run.

Checking Equilibriation

The simulation takes about 0.3 seconds to equilibriate energy, temperature, and pressure, as shown below:

Energy equilbriation

Closer look at energy equilibriation

Temperature equilbriation

Closer look at temperature equilibriation

Pressure equilbriation

Closer look at pressure equilibriation

Energy equilibriation for multiple timesteps

Of the five timsteps used, 0.0025 is the largest acceptable timestep to use as a smaller timestep of 0.001 results in a very similar equilibriation, so going this small is not necessary. 0.015 does not equilbriate at all as the time steps are too large for the numerical integration to accurately find an average for the ensemble, and energy drifts; diverging instead of converging to an average value.

Running Simulations Under Specific Conditions

Thermostats & Barostats

$E_K=\frac{3}{2}N{k}_B\mathfrak{𝔗}=\frac{1}{2}\sum _i{m}_i(\gamma v_i{)}^2$

${\gamma }^2=\frac{\frac{3}{2}N{k}_B\mathfrak{𝔗}}{\frac{1}{2}\sum _i{m}_i{v}_i^2}$

${\gamma }^2=\frac{\frac{3}{2}N{k}_B\mathfrak{𝔗}}{\frac{3}{2}N{k}_{B}T}$

$\gamma =(\frac{\mathfrak{𝔗}}{T}{)}^{\frac{1}{2}}$

Examining the Input Script

The command fix aves all ave/time 100 1000 100000 means that values will be sampled every 100 timesteps; in total 1000 readings will be taken to compute a final average on the 100000th timestep.

run 100000 indicates that 100000 timesteps will be simulated.

Plotting the Equations of State

A plot of density against temperature compared to ideal gas law

Higher pressures lead to higher densities, both in theory and in these simulations. Our simulated density is higher than that given by the ideal gas law because the simulation takes particle interactions into account. The error increases at higher pressures, when more collisions are likely to occur, while lower pressures would theoretically behave more as an ideal gas would. For the same reason, error decreases at higher temperatures.

Calculating Heat Capacities Using Statistical Physics

A plot of heat capacity/volume against temperature

Higher pressure results in higher heat capacity as the increased number of molecules per unit volume that can absorb energy to their vibrational excited states. As the simulation is in a lattice, rotational degrees of freedom are not available to the atoms, and so heat capacity decreases as temperature increases, despite expectations. Alternatively, increasing temperature causes a decrease in band gap, requiring less energy to enter excited states and thus lowering heat capacity.

An example of the input scripts is below:

Media:inputnpt.in

Structural Properties and the Radial Distribution Function

A plot of RDFs

The RDF shows the probability of finding a particle at a distance r from a reference particle, relative to an ideal gas. In a gas, there is little order and minimal structure to particles and so the graph has minimal features.

Liquids are slightly more ordered and the decreasing heights of peaks of the RDF correlate to coordination spheres. There is a high probability of finding another particle in a primary coordination sphere but this probability decreases as you go farther away from the reference particle, indicated by decreasing heights of peaks.

The solid FCC lattice has a much higher order, and the RDF peak separation and heights define the lattice structure. The first, second, and third sharp peaks refer to different sets of nearest neighbours, while their heights show how many of those nearest neighbours there are. The lattice spacing is the same as the distance to the second nearest neighbour, 1.475. This agrees well with the original input density of 1.3 (which should result in a lattice spacing of 1.45).

A diagram of first (left) and second (right) nearest neighbours (shown in blue) with respect to a reference particle (red) in an FCC lattice (other atoms shown in black)

A diagram of third nearest neighbours (blue) with respect to a reference particle (red) in an FCC lattice (other atoms shown in black)

A plot of the running integral of the RDF of a solid

The coordination numbers are 12 (Int(g(1.205)=12, 12 neighbours), 6 (Int(g(1.475)=18, 18-12=6 neighbours), and 24 (Int(g(1.775)=42, 42-18=24 neighbours) respectively.

Dynamical Properties and the Diffusion Coefficient

Mean Squared Displacement

• $D_{solid}=\frac{1}{6}\cdot 0.00382752458475=6.38x10^{-4}$
• $D_{liquid}=\frac{1}{6}\cdot 0.509774821123=0.085$
• $D_{gas}=\frac{1}{6}\cdot 7.99193006423=1.33$

The diffusion coefficient increases as entropy of the phase increases, which matches expectations as gas particles are much more likely to diffuse than a rigid lattice of solid molecules.

Total MSD as a function of time of simulated solid

Total MSD as a function of time of simulated liquid

Total MSD as a function of time of simulated gas

1000000 Atoms

Total MSD as a function of time for 1000000 atoms

• $D_{solid}=\frac{1}{6}\cdot 2.79195534196x10^5=4.65x10^{-6}$
• $D_{liquid}=\frac{1}{6}\cdot 0.531614512766=0.0886$
• $D_{gas}=\frac{1}{6}\cdot 18.0968139669=3.02$

The MSD graph for a gas is curved at first, indicating ballistic motion proportional to $T^2$. After enough collisions have occurred, diffusion is linear, as it is for a liquid which constantly has the same collisions. The diffusion coefficient is close to 0 for solids which is as expected.

Velocity Autocorrelation Function

• $x(t)=Acos(\omega t+\varphi )$ and $v(t)=\frac{dx(t)}{dt}=-A\omega sin(\omega t+\varphi )$

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}v(t)-v(t+\tau )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{v}^2(t)dt}=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )sin(\omega (t+\tau )+\varphi )}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )dt}$

• $sin(A+B)=sinAcosB+cosAsinB$

$C(\tau )=\frac{[cos(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )dt]+[sin(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )cos(\omega t+\varphi )dt]}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )dt}$

• $si{n}^2(x)=\frac{1}{2}(1-cos(2x))$

$C(\tau )=\frac{[cos(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{1}{2}(1-cos(2(\omega t+\varphi )]+[sin(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )cos(\omega t+\varphi )dt]}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{1}{2}(1-cos(2(\omega t+\varphi )dt}$

• $\int \frac{1}{2}(1-cos(2(\omega t+\varphi ))dt=\frac{t}{2}+\frac{sin(2(\omega t+\varphi )}{4\omega }+c$
• $\int sin(\omega t+\varphi )cos(\omega t+\varphi )dt=\frac{si{n}^2(\omega t+\varphi )}{2\omega }+c$

$C(\tau )=\frac{{[\frac{tcos(\omega \tau )}{2}+\frac{cos(\omega \tau )sin(2(\omega t+\varphi ))}{4\omega }+\frac{sin(\omega \tau )si{n}^2(\omega t+\varphi )}{2\omega }]}_{-\mathrm{\infty }}^{\mathrm{\infty }}}{{[\frac{t}{2}+\frac{sin(2(\omega t+\varphi )}{4\omega }]}_{-\mathrm{\infty }}^{\mathrm{\infty }}}$

• sin(x) is an odd function and integrating between $-\mathrm{\infty }$ and $\mathrm{\infty }$ will result in 0

$C(\tau )=\frac{\frac{tcos(\omega \tau )}{2}}{\frac{t}{2}}=cos(\omega \tau )$

VACF minima refer to collisions of particles where velocity is instantaneously 0, negative as they are in the opposite direction. As VACF is averaged over all molecules, they cancel out once they are out of phase, which happens faster for liquids than it does for solids. In comparison to the harmonic oscillator, which only models one particle without any collisions, no convergence to 0 occurs.

A plot of total VACF against $\tau$

Diffusion coefficient estimations, using the trapezium rule:

• $D_{solid}=\frac{1}{3}\cdot 0.43184744815700105=0.144$
• $D_{liquid}=\frac{1}{3}\cdot 146.83331703729999=48.9$
• $D_{gas}=\frac{1}{3}\cdot 1451.848385=484$

Plot of running integral of VACF of a solid

Plot of running integral of VACF of a liquid

Plot of running integral of VACF of a gas

1000000 Atoms

A plot of total VACF against $\tau$ for 1000000 atoms

• The trapezium rule estimation of the integral for a solid was found to be -0.416 for 1000000 atoms between 0 and 500. (D would hypothetically equal -0.139)
• $D_{liquid}=\frac{1}{3}\cdot 123.7270701106=41.2$
• $D_{gas}=\frac{1}{3}\cdot 1466.443215=489$

Plot of running integral of VACF of a solid

Plot of running integral of VACF of a solid

Plot of running integral of VACF of a solid

The estimated diffusion coefficients for the two simulations follow the same trend, however the values obtained for the larger 1000000 atom simulations are generally smaller. The largest sources of error include the trapezium rules used to calculate the integral and the simulation assumption that velocities do not change upon collisions.