Rep:Mod:ZG Liquid Simulation

Introduction to Molecular Dynamics Calculations

Harmonic Oscillator

Calculating the difference in position between the Velocity-Verlet (V-V) algorithm and analytical solution, seen in equation 1 and 2 respectively (where $x$, $v$,$t$, $\delta t$, $A$, $\omega$, $\varphi$ are position, velocity, time, timestep, amplitude of vibration, frequency and phase shift.), revealed the accuracy of the V-V algorithm. Taking the modulus of this difference ensured that no negative values were obtained.

$x_i(t+\delta t)=x_i(t)+v_i(t+\frac{1}{2}\delta t)$ (1)

$x(t)=Acos(\omega t+\varphi )$ (2)

The energy of a classical simple harmonic oscillators described by equation 3; where $E$, $U$, $K$, $m$, $v$ and $k$ are the total energy, potential energy, kinetic energy, mass, velocity and spring constant, respectively.

$E=U+K=\frac{1}{2}m{v}^2+\frac{1}{2}k{x}^2=\frac{1}{2}k{A}^2$ (3)

Substitution of position and velocity, from the V-V algorithm, seen in equation 1 and 4 (where $a$ is acceleration) respectively, into equation 3 permits the total energy of the harmonic oscillator to be determined from this approximate method. As the oscillator exhibits simple harmonic motion, there is no damping forces to remove energy from the system and therefore, no reduction in the amplitude occurs; it is assumed that the spring constant of the oscillator does not change either. Thus, total energy should remain constant.

$v_i(t+\delta t)=v_i(t+\frac{1}{2}\delta t)+\frac{1}{2}{a}_i(t+\delta t)$ (4)

Error Plot

Figure 1 displays the difference between calculated (from the V-V theorem with a timestep of 0.1) and analytical solution as a function of time; this difference has been referred to as the error. It is evident that periodic oscillations in the error occurred, with increasing maximums in the error as time increased.

Figure 1. Error as a function of time for a timestep of 0.1 s between analytical and Velocity-Verlet algorithm.

If the derivative of the error is taken with respect to time, the time of maximum errors can be extracted from the points closest to zero. Figure 2 displays the maximum error as a function of time. It is evident that there is a linear increase in maximum error with time, which shows that the uncertainty in the V-V algorithm from analytical increases with time.

Figure 2. Fitted maximum error as a function of time for a timestep of 0.1 s between analytical and Velocity-Verlet algorithm.

Every calculation performed had a small, but finite error, which accumulates to significant values with time. An equation can be obtained for the linear fit which allowed the error for each calculation and errors for specified time periods to be calculated. The error in each calculation was found to be $4\times 10^{-5}$, which was obtained from multiplying the timestep by the gradient of the linear function. A value of the error associated to a time period, which was taken to be 1 s for simplicity, was obtained by multiplying the error in each calculation by the number of calculations in the specified time period; for 1 s, this had a value of $4\times 10^{-4}$. A time period of 2500 s would have to pass before a 1 % error arose.

Timestep

A time step of 0.2 s produced an error of 1% in the total energy between analytical and V-V algorithm. Monitoring total energy of the system ensures no divergences from analytical. For a closed system, such as a simple harmonic oscillator, the total energy should remain constant, as seen by equation 3. Clearly, if the system progressively loses energy, then erroneous calculations are being performed. As the V-V algorithm is an approximation, it would not be expected that the energy exactly matched that of the analytical solution. With an increased timestep, the approximation becomes larger, which causes larger changes in the energy.

Figure 3. Energy as a function of time for the Velocity-Verlet algorithm with a timestep of 0.2 s.

A plot of the total energy, for a time step of 0.2 s, can be seen in Figure 3. Periodic changes can clearly be seen in the total energy, which do not occur for the analytical solution.

Lennard-Jones

The separation at zero potential energy $(r_0)$ has a value of $\sigma$. At infinite separations the potential approaches zero as well, but this value has no importance.

Force (F) can be calculated from the negative derivative of potential with respect to position. This calculation was performed on the Lennard-Jones (L-J) potential, which is displayed in equation 5, and presented in equation 6; where $\varphi$, $r$, $ϵ$ and $\sigma$ are the potential energy, separation, depth of potential well and interatomic distance at which potential reaches zero, respectively. Calculation of the force ($F(r_0)$) at $(r_0)$ can be found by substitution the result from earlier; the force at zero potential is displayed in equation 7.

$\varphi (r)=4ϵ\left(\right(\frac{\sigma }{r}{)}^{12}-\left(\frac{\sigma }{r}{)}^6\right)$ (5)

$F(r)=-\frac{d\varphi (r)}{dr}=\frac{24ϵ}{\sigma }\left(2\right(\frac{\sigma }{r}{)}^{13}-\left(\frac{\sigma }{r}{)}^7\right)$ (6)

$F(r_0)=\frac{24ϵ}{\sigma }$ (7)

Equilibrium separation occurs when the interatomic force acting equals zero, which occurs at the minimum of the L-J potential. Equation 6 can be equated to zero and solved for $r$ to yield $r_{eq}$; the result for which is displayed in equation 7. Well depth was found by substituting $r_{eq}$ into the L-J potential, which was found to have a value of $ϵ$. Hence, to break an atomic bond, in the L-J potential, $ϵ$ must be supplied.

$F(r)=0=\frac{24ϵ}{\sigma }\left(2\right(\frac{\sigma }{r}{)}^{13}-\left(\frac{\sigma }{r}{)}^7\right)$

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

$r_{eq}=\sqrt[6]{2}\sigma$ (8)

An indefinite integral was taken of the L-J potential, so to allow substitution of limits for each question. The result of the indefinite integral can be seen in equation 9.

$\int \varphi (r)dr=4ϵ\sigma \left(\frac{1}{5}\right(\frac{\sigma }{r}{)}^5-\frac{1}{11}\left(\frac{\sigma }{r}{)}^{11}\right)+c$ (9)

Each integral can be solved by simply substituting the relevant limits into equation 9, which has been carried out and displayed in equations 10 through 12.

$\underset{\mathrm{\infty }}{\overset{2\sigma }{\int }}\varphi (r)dr=4ϵ\sigma \left(\frac{1}{5}\right(\frac{\sigma }{r}{)}^5-\frac{1}{11}\left(\frac{\sigma }{r}{)}^{11}\right){|}_{\mathrm{\infty }}^{2\sigma }=\sigma ϵ(0-0.0248)=-2.48\times 10^{-02}$ (10)

$\underset{\mathrm{\infty }}{\overset{2.5\sigma }{\int }}\varphi (r)dr=4ϵ\sigma \left(\frac{1}{5}\right(\frac{\sigma }{r}{)}^5-\frac{1}{11}\left(\frac{\sigma }{r}{)}^{11}\right){|}_{\mathrm{\infty }}^{2.5\sigma }=\sigma ϵ(0-8.18\times 10^{-03})=-8.18\times 10^{-03}$ (11)

$\underset{\mathrm{\infty }}{\overset{3\sigma }{\int }}\varphi (r)dr=4ϵ\sigma \left(\frac{1}{5}\right(\frac{\sigma }{r}{)}^5-\frac{1}{11}\left(\frac{\sigma }{r}{)}^{11}\right){|}_{\mathrm{\infty }}^{3\sigma }=\sigma ϵ(0-3.29\times 10^{-03})=-3.29\times 10^{-03}$ (12)

From inspection of the L-J potential, it would be expected that the area obtained is negative and decreases with increasing $\sigma$, as the potential is negative and approaches zero at large interatomic distances.

It can be seen that the area from 2 sigma to infinity is an order of magnitude larger than that of 3 sigma to infinity. Hence, under the error of the experiment, it would be safe to perform calculations up to, or a little further than, 3 sigma instead of to infinity. Thereby, reducing the computational expenditure.

Atomic Number Estimation

The molar volume of water, under standard conditions, is 18 ml. Dividing 1 ml by the molar volume yields the number of moles of water in the specified volume; multiplying by Avogadro’s number ($6.022{\times }^{23}$) affords the number of atoms, which turned out to be $3.34{\times }^{22}$. Conversely, dividing the number of atoms by Avogadro’s number and multiplying by the molar volume of water produces its volume, which for 1,000 atoms of water produced a volume of $2.99{\times }^{-20}$ ml. Clearly, there is significant differences between the volumes and number of atoms in each case. As the number of atoms increases in a molecular dynamics calculation, so does the time taken, owing to a larger number of calculations being performed. Hence, relatively small numbers of atoms are taken for calculations, instead of what is considered a relatively small volume (1 ml).

Periodic Boundary Conditions

After applying periodic boundary conditions, the atom would have coordinates (0.2, 0.1, 0.7). The atom reaches the limit of the x-value first, passes through the yz-plane and reappears at (0, 13/14, 9/14); then the atom reaches the limit of the y-value, and passes through the xz-plane at (7/84, 1, 56/84) and reappears at (7/84, 0, 56/84).

Reduced Units

The Lennard-Jones cutoff is r = 10.88 nm in SI units. In SI units, the well depth is - 0.998 KJ mol^-1. A reduced temperature of T* = 1.5 has a value of 180 K in SI units.

Equilibration

Random Coordinates

In a liquid, there is no periodic long range order, so it would not be expected that all interatomic distances are equal. Hence, random generation is a possibility to produce atomic coordinates. However, there is one main problem with generating random coordinates for atoms in a liquid; that being: atomic coordinates that do not correspond to realistic interatomic distances, be it too small or large. If coordinates are generated randomly, it is possible for two atoms to have exactly the same position. This is the extreme case, and is unlikely to occur, but in no circumstance could this occur under normal laboratory conditions. More likely, and less extreme, is the generation of coordinates that are smaller than typical interatomic distances. Therefore, increasing the energy of the system beyond what would be expected. On the other hand, regions of space could exist where atoms do not reside at all, but might be expected to. This is the main reason not using randomly generated atomic coordinates.

Lattice points

A primitive cubic lattice contains one lattice point in the unit cell. The volume of the unit cell, for a cubic system, can simply be calculated by cubing one of the lattice parameters. Dividing the number of lattice points by the primitive unit cell volume yields the number density. If this is done for a lattice parameter of 1.07722, then a number density of 0.8 is obtained. A face-centered cubic unit cell contains 4 lattice points; the contribution from corner atoms equals one and those from faces three. The volume of the unit cell can be obtained by dividing the number of lattice points by the density of points; taking the cube root of this value yields that lattice parameter of the unit cell, which had a value of 1.49.

Atoms from Lattice

Face-centered cubic lattices have four lattice points on the unit cell; one from atoms positioned at corners and three from atoms residing on faces. As the box contains 1000 unit cells of the lattice there would be 4000 atoms.

LAMMPS manual

The first lines specifies the mass of atoms. In this line, one type of atom exists with a numerical value of 1. Specification of the cutoff point of pairwise interactions, for the L-J potential, is provided in the second line, which has been set equal to $3\sigma$. This was determined in a previous question to be a distance at which no significant interaction occurred. Finally, specification of pairwise force field coefficients is performed. A leading and trailing asterisk indicates that all atoms should have the following parameters. The L-J potential requires two force field constants: $\sigma$ and $ϵ$; both have been set to one.

Integration Algorithm

Velocity-Verlet algorithm, the initial conditions for which are $x_i(0)$, $v_i(0)$, as specified.

Running Simulation

These two lines are the final statement lines of each section of code. No information is provided about what calculation is being performed, or with what parameters. Despite the lack of information, these lines should still proceed with the same calculation.

Timestep has not been verified as a variable, so it can not be used in a subsequent calculation. Hence, the code would not know what timestep to take when performing the V-V algorithm, if it was prompted with a timestep variable. Furthermore, the number of steps was not determined from dividing calculation duration by timestep. Therefore, the number of steps can not be used in proceeding parts of the code and the duration of the calculation has not been explicitly stated.

Checking Equilibration

Plots of the converged reduced energy, reduced temperature and pressure, as a function of time, can be seen in figures 4 through 6, respectively. In all cases, for the 0.001 timestep, it was found that the calculations converged. Reduced energy, reduced temperature and pressure converged to a mean value of -3.18, 1.25, and 2.60, respectively. From inspection of the figures, it is evident that convergence of the parameters was reached by 0.4 s. Random fluctuations occurred after this time, but this was expected and is not to be of concern.

Figure 4. Reduced energy as a function of time for a timestep of 0.001 s over a 1 s time period.

Figure 5. Reduced temperature as a function of time for a timestep of 0.001 s over a 1 s time period.

Figure 6. Pressure as a function of time for a timestep of 0.001 s over a 1 s time period.

Figure 4 displays reduced energy converging to its mean value without overshooting and with the smallest fluctuations out of these parameters. Random variations in reduced temperature and pressure are more significant and both initially overshot their equilibrium value. The reduced energy is already close to its converged value, but reduced temperature and pressure are not. Hence, changes in reduced temperature and pressure are initially large and their values overshoot the equilibrium before converging.

Comparing figure 5 and 6 revealed a link between reduced temperature and pressure. Initially the pressure increases, but the reduced temperature decreases. Then the pressure overshoots at more positive values than its equilibrium and subsides to its converged equilibrium; whereas, reduced temperature decreases to more negative values than equilibrium. Hence, it can be inferred, from figures 5 and 6, that the pressure and temperature are inversely proportional, which would be expected for other theories, such as the ideal gas law. If fluctuations at times larger than 0.4 s are analysed, it can be seen that the trend continues.

Figure 7 displays energy as a function of time for all of the calculated timesteps. Clearly, a timestep of 0.015 does not converge to a value, and therefore, can not be used for any calculation. Timesteps of 0.01 and 0.0075 converged, but to values that were slightly above that produced by a timestep of 0.001. They could be used as timesteps, but would produce results that were slightly erroneous. Both 0.001 and 0.0025 timesteps essentially converged to the same energy value. As the latter requires less computational time, it shall be used as a compromise between accuracy and efficiency.

Figure 7. Reduced energy as a function of time for timesteps 0.001, 0.0025, 0.0075, 0.01 and 0.015, which can be seen in light blue, orange, grey, yellow and dark blue, for a time period of 100 s.

Running Simulations Under Specific Conditions

Thermostats and Barostats

Owing to the fact that $\gamma$ is not part of the summation, it can be pulled out and solved for. Both equations were solved for $3N{k}_b$, seen in equations 13 and 14, and equated; the resulting equation was solved for $\gamma$, which yielded equation 15.

$\frac{\sum _i{m}_i{v}_i^2}{T}=3N{k}_b$ (13)

$\frac{{\gamma }^2\sum _i{m}_i{v}_i^2}{\tau }=3N{k}_b$ (14)

$\frac{{\gamma }^2\sum _i{m}_i{v}_i^2}{\tau }=\frac{\sum _i{m}_i{v}_i^2}{T}$

$\gamma =\sqrt{\frac{\tau }{T}}$ (15)

If equation 15 is inspected, it can be seen to hold true for its dependence on the relative magnitudes of $T$ and $\tau$. For $T>\tau$ it would be expected that $\gamma$ is smaller than 1.

Importance of the Three Numbers

The first number its referred to as Nevery, which determines how frequently an input value should be utilised for the average. It has a value of 100 in this example; hence, every 100th input value is taken and used in the subsequent calculations. How frequently values were taken, in terms of time, can be determined by multiplying Nevery with the timestep. Hence, values are taken every 0.1 s, as there is a 0.001 s timestep and 100 steps between values being used.

The second number, with a value of 1,000, referred to as Nrepeat, specifies the number of input values were used for determining averages. A total time until an average is determined can be found from multiplying the frequency of sampling, which was determined in the previous paragraph, by the number of values used to determine the average. Hence, the product of timestep, Never and Nrepeat yields the total time, which had a value of 100 s for this example.

Finally, Nfreq sets the number values timesteps that pass before the average is calculated. Nfreq can be found from taking the product of Nevery and Nrepeat.

Equations of State

Figure 8 displays density as a function of reduced temperature for two pressures from molecular dynamics calculations and predicted values from the ideal gas law. Error bars were obtained from the standard deviations in density and temperature, which were displayed in the .log file.

It is evident that the ideal gas law predicts significantly higher density values than the V-V algorithm. If the assumptions of the ideal gas model are considered, then this observation can be rationalised. The ideal gas model assumes that all atoms can be treated as point like particles that do not interact. Hence, when pressure is applied, no repulsive forces arise and atoms approach without a lower limit. Whereas, molecular dynamics models interactions of atoms from interatomic potentials, such as the L-J potential, which was used here. Once an internuclear separation of $\sigma$ is reached, atoms start to repeal. As the interatomic distances reduces further, the potential increases proportionally to $1/r^{12}$, which ensures that particles do not reach internuclear distances significantly shorter than $\sigma$. A more accurate description of the density could be obtained with the Van der Waals equation, owing to the fact that this model takes into account the atomic interactions and volumes.

Figure 8. Density as a function of reduced temperature for two pressures, 2.0 and 3.0 seen in blue and red respectively, with Velocity-Verlet algorithm calculations and the predicted density from the ideal gas law, seen in green and yellow, for pressures of 2.0 and 3.0, respectively.

As reduced temperature increases the difference between calculation methods reduces, which occurred because of the inverse relationship for the ideal gas law initially reducing at a larger rate than the linear V-V algorithm.

It can be seen that a larger pressure results in a higher density, which is intuitively expected for most materials.

Heat Capacity

Heat capacity over volume was calculated from the adjustments made to the script, which is displayed below.

### DEFINE SIMULATION BOX GEOMETRY ###
variable dens equal 0.8 ### This line was added so the density could be explicitly changed 
lattice sc ${dens} ### The density of the lattice is specified here 
region box block 0 15 0 15 0 15
create_box 1 box
create_atoms 1 box

### DEFINE PHYSICAL PROPERTIES OF ATOMS ###
mass 1 1.0
pair_style lj/cut/opt 3.0
pair_coeff 1 1 1.0 1.0
neighbor 2.0 bin

### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###
variable T equal 2.0  ### The temperature was varied here
variable timestep equal 0.0025 ### A timestep of 0.0025 s was found to be sufficient 

### ASSIGN ATOMIC VELOCITIES ###
velocity all create ${T} 12345 dist gaussian rot yes mom yes

### SPECIFY ENSEMBLE ###
timestep ${timestep}
fix nve all nve

### THERMODYNAMIC OUTPUT CONTROL ###
thermo_style custom time etotal temp
thermo 10

### RECORD TRAJECTORY ###
dump traj all custom 1000 output-1 id x y z

### SPECIFY TIMESTEP ###

### RUN SIMULATION TO MELT CRYSTAL ###
run 10000
unfix nve
reset_timestep 0

### BRING SYSTEM TO REQUIRED STATE ###
variable tdamp equal ${timestep}*100
fix nvt all nvt temp ${T} ${T} ${tdamp} ### Changed as described 
run 10000
reset_timestep 0

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp density
variable erg equal etotal
variable erg2 equal etotal*etotal
variable temp equal temp
fix aves all ave/time 100 1000 100000 v_temp v_erg v_erg2
### Pressure and density calculations were removed for efficiency 
unfix nvt ### Changed as described 
fix nve all nve ### Changed as described 
run 100000  

variable avetemp equal f_aves[1]
variable heatc equal (atoms*atoms*(f_aves[3]-f_aves[2]*f_aves[2]))/(${avetemp}*${avetemp}*vol)
### The line above calculates the heat capacity per unit volume 

print "Averages"
print "--------"
print "Temperature: ${avetemp}"
print "Heat Capacity over Volume: ${heatc}"

The heat capacity over volume was calculated for two density values and five reduced temperatures; figure 9 displays the obtained results.

Figure 9. Heat capacity over volume as a function of reduced temperature for densities of 0.2 and 0.8, seen in blue and red, respectively.

Figure 10. Heat capacity over volume as a function of the reciprocal of the reduced temperature cubed for densities of 0.2 and 0.8, seen in blue and red, respectively.

Heat capacity is defined as the energy required per unit increase in temperature, for constant volume in this case. A reduction in heat capacity would be expected for a solid, at low temperatures, after the Schottky anomaly. Heat capacities tend to increase, or stay the same, with respect to temperature, owing to the fact that, at higher energies, the difference in energy between states is lower; hence, less energy is required for excitation, which causes larger changes in temperature per unit energy supplied. There comes a point however, when atoms posses enough energy to interchange between states without the need for extra energy put, owing to the fact that the energy levels do not have large differences and atoms can interchange energy. This causes the system to not absorb energy being supplied, so the heat capacity starts to decrease.

Thermal expansion of phases tends to occur with increasing temperature. Owing to the fact that the heat capacity has been divided by volume, the reduction in heat capacity over volume, with respect to temperature, was exaggerated. Expansions occur at larger temperatures because atoms have larger translational and vibrational motion.

It would be expected that a higher density corresponds to a larger heat capacity over volume, owing to the larger number of atoms per unit volume.

Clearly, both densities have lower heat capacities over volume for larger reduced temperatures. Equation 16 displays the relationship between heat capacity over volume and temperature, where the ideal gas law has been substituted to include the density dependence. This equation revealed that the heat capacity over volume is directly proportional to the reciprocal of the reduced temperature cubed. Hence, plotting the heat capacity over volume against the reciprocal of the reduced temperature cubed should yield a linear relationship, which can be seen in figure 10.

$\frac{C_v}{V}=Var[E]\frac{NP}{k_b^2{T}^3}$ (16)

Radial Distribution Function

Figure 11 displays the radial distribution function (RDF) as a function of reduced radius for the solid, liquid and vapour phase.

Figure 11. G(r) as a function of reduced radius for the solid, liquid and vapour phase, which are seen in blue, orange and green, respectively.

Figure 12. Running integral of g(r) as a function of reduced radius for the solid, liquid and vapour phase, which are seen in blue, orange and green, respectively.

All phases have zero $g(r)$ up to one, which is expected, as these distances correspond to effective interatomic bond distances and therefore, no atoms should reside here. Furthermore, all phases increase suddenly, from this zero point, and reach a maximum before subsiding to value near, or oscillating about, one. The solid phase has the largest maximum point in $g(r)$, followed by the liquid phase and vapour phase, which has the smallest maximum. These differences in $g(r)$ would be expected between the phases, as the solid phase is slightly denser than liquid, both of which are significantly denser than the vapour phase. All phases converge to one as their RDF values have been normalised to the bulk value.

Both solid and liquid phases exhibit oscillatory changes in $g(r)$ after reaching a maximum shortly after $\sigma$; this behaviour is not observed for the vapour phase. Oscillatory nature in the RDF indicates that these phases pack in such a way that they form “shells” and that these atoms tend to remain in constant positions relative to adjacent atoms. As there is no order in the vapour phase, the formation of “shells” around atoms is not seen. Furthermore, atoms in the vapour phase constantly move relative to each other, which prevents spikes in the RDF. Larger oscillatory changes in the RDF of the solid phase, when compared to the liquid phase, are observed because of regularly packed stationary atoms. Whereas, in liquid phase, atoms are constantly moving, so dynamic “shells” of atoms form.

Interatomic distance at equilibrium for atoms described by the L-J potential is given by $\sqrt[6]{2}$ (1.12), in reduced units. There is only one maximum in the vapour RDF curve, which occurs at a reduced radius of 1.125. Initial, maximums in the RDF of the liquid and solid phase occurred at 1.25 and 1.08, respectively; hence, it can be concluded that all of these values are in relatively good agreement with the expected.

Subsequent maximums, in the solid phase, occurred at reduces radii of 1.58 and 1.98. The first maximum in the RDF, with a reduced radius of 1.25, corresponds to the distance between a corner atom and an atom located on a face, which has a corner on the corner atom. Atoms on adjacent corners produce the next closest maximum (1.58). An interatomic maximum is produced at a reduced radius of 1.98 from the distance between an atom on a corner and an atom on the face of the unit cell that does not share a corner with the corner atom. Figure 13 displays these relationships.

Figure 13. Face-centered cubic unti cell decorated with the three maximums that occur at the shortest reduced radius in the RDF.

Assuming that the shortest distance is between a corner atom and one located on the face that shares a vertex with the corner atom, the distances to all other atoms can be inferred. All of these calculations were completed with simple trigonometry. Equation 17 mathematically describes the relationship between A and B, from figure 13.

$B=2ACos(45^o)=\frac{2A}{\sqrt{2}}$ (17)

Substituting the values yields a value of 1.59 for B, which is in good agreement with the obtained peek from the RDF. Calculation of the next peek can be achieved from pythagorus' theorem, seen in equation 18.

$C=\sqrt{B^2+(B/2{)}^2+(B/2{)}^2}=\sqrt{\frac{3}{2}}B=\sqrt{3}A$ (18)

Again the calculated is in good agreement the obtained for the RDF, with a value of 1.95 being obtained for C.

A running integral of $g(r)$ was calculated, and has been displayed in figure 12. It would be expected that the solid phase has the largest area underneath the $g(r)$, owing to the fact that it is the phase with the largest density. This expectation is seen in figure 12, which also accurately predicts that the liquid has the next largest density. Integration of peaks, in the solid phase, revealed the number of atoms which hold the relationships described earlier. Taking the range of reduced radii between 0.925 and 1.375, the integration comes to a value of just over 12. This would be expected, as it is known, for a face-centered cubic lattice type, that each atom has 12 nearest neighbours. Integrating the peak between reduced radii of 1.375 and 1.725 yields a result of 6, which indicates that there are 6 atoms of a B distance. Finally, integrating from 1.825 through 2.125 produces a area of 24. Hence, there are 24 atoms that reside a distance of C from each atom.

As reduced radius increases, the running integral should become linear as the density approaches the bulk value.

Dynamical Properties

From the phase diagram for the L-J potential, it was found that a temperate of 2.0 and density of 1.5 should produce a solid phase; the vapour phase was produced by a reduced temperature of 2.0 and density of 2.0. Initially, a lower reduced temperature (0.5) was used to estimate the solid phase. However, it was found that the output files resembled the liquid phase more than solid, when compared to calculations with 1,000,000 atoms.

Mean Square Displacement

Owing to the fact that the order of atoms decreases as the phases changes form solid to a vapour, it would be expected that the mean square displacement of a vapour would be larger than a solid. Figure 14 displays the mean square radius as a function of timestep for 8,000 atoms. It is evident, from this figure, that the vapour phase, seen in green, has a significantly larger mean square radius than any other phase. Owing to the fact that atoms in the vapour phase have large translational velocities and weak interatomic interactions, this observation would be expected. As a liquid has a similar density to a solid, but no long range order that anchors atoms in fixed locations, it would be expected that the liquid phase has a larger mean square displacement. This expectation is seen in figure 14 by the orange line reaching higher values of the mean square displacement than blue.

Figure 14. Mean square displacement as a function of timestep for solid, liquid and vapour phases, seen in blue, orange and green, respectively.

Figure 15. Diffusion coefficient as a function of time for solid, liquid and vapour phases, seen in blue, orange and green, respectively.

Figure 16 also displays the mean square displacement as a function of timestep, but for one million atoms. If figure 14 and 16 are compared, it is evident that they are essentially identical. The same observations for each phase are consistent. Hence, it can be concluded that 8,000 atoms is sufficient to perform these calculations, and should be used preferentially over 1,000,000 atoms as it is less computationally expensive.

Equation 19 was utilised to calculate to diffusion coefficient of the atoms in each phase. Figures 15 and 17 display the diffusion coefficient as a function of time for 8,000 and 1,000,000 atoms, respectively. Again, there is essentially no difference between the graphs.

It would be expected, just like in the mean square displacement, that the vapour phase has the largest diffusion coefficient and the solid phase has the smallest; this is, in fact, seen in these figures. Atoms in the vapour phase have significant translational velocities and weak interatomic forces compared to the other phases. Furthermore, spaces between atoms, in the vapour phase, is orders of magnitude larger than the solid or liquid. Hence, the atoms can translate with being impeded by other atoms.

$D=\frac{1}{6}\frac{\partial <r^2(t)>}{\partial r}$ (19)

Figure 16. Mean square displacement as a function of timestep for solid, liquid and vapour phases, seen in blue, orange and green, respectively. These values were obtained for the provide files with 1,000,000 atoms.

Figure 17. Diffusion coefficient as a function of time for solid, liquid and vapour phases, seen in blue, orange and green, respectively. These values were obtained for the provide files with 1,000,000 atoms.

The diffusion coefficient in the vapour phase increases over time and converges to a value of 3.20 for 8,000 atoms. In a previous question, it was noted that the initial positions of atoms are obtained from positioning atoms on lattice points of a simple cubic lattice type. Hence, the density of the phase is larger, which to some extent, prevents atoms from moving large distances without colliding and interacting with other atoms. Furthermore, atoms initially have zero velocity, from the starting condition of the V-V algorithm, which means their displacements are limited.

The table below displays the calculated diffusion coefficients for 8,000 and 1,000,000 atoms in the solid, liquid and vapour phase. There is very little difference between the calculated diffusion coefficient, for each phase, with different numbers of atoms. As expected, and seen in figures 15 and 17, there are large differences between the phases.

Phase	D / cm2s-1	D / cm2s-1 (Million)
Solid	$8.30\times 10^{-06}$	$2.50\times 10^{-05}$
Liquid	$8.90\times 10^{-02}$	$8.90\times 10^{-02}$
Vapour	$3.20$	$3.30$

Velocity Autocorrelation Function

The analytical solution of position for the simple harmonic oscillator can be seen in equation 2. A derivative of position was taken with respect to time to yield velocity of the harmonic oscillator as a function of time; equation 20 displays the solution.

$v(t)=\frac{dx(t)}{dt}=-A\omega sin(\omega t+\varphi )$ (20)

Substitution equation 20 into the integral for the normalised velocity autocorrect function for a 1D harmonic oscillator yields equation 21.

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

The addition formula of sin can be utilised for the $t+\tau$ sin function to simplify the integral. Equation 22 displays this substitution.

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

Simplification of equation 22 can be achieved by splitting the integral into two separate fractions, which is displayed in equation 23.

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

Owing to the fact that the integral of an odd function, between two limits that are equidistant from the origin, is equal to zero, the simplification furthers to equation 24.

$C(\tau )=cos(\omega \tau )$ (24)

Figures 18 and 19 display the VACF calculations as a function of time for 8,000 and 1,000,000 atoms in the liquid and solid respectively, with superimposed analytical solution for the 1D harmonic oscillator.

Figure 18. VACF as a function of time for solid and liquid phases, seen in blue and orange respectively, plotted with the analytical solution for the 1D HO, which is seen in grey.

Figure 19. VACF as a function of time for solid and liquid phases, seen in blue and orange respectively, plotted with the analytical solution for the 1D HO, which is seen in grey. These values were obtained for the provide files with 1,000,000 atoms.

There is essentially no minima for the liquid phase. However, the solid phase has one major minimum point, which occurred shortly after the initial conditions, and several smaller minimums. For 8,000 atoms, it can be seen that the curative of the minimum resembles that of the 1D harmonic oscillator. Hence, it could be suggested that is simple harmonic motion occurred. After this initial minimum, the frequency of these vibrations, for 8,000 atoms, is very similar to that of the 1D HO. However, for one million atoms, the frequency of these vibrations is roughly twice that of the 1D HO.

It can be seen in figure 18 and 19 that there is a difference between the solid and liquid phase near initial conditions; as time increases, both phases converge to zero. For the liquid phase, there is no periodic changes in the velocity autocorrelation function, it exponentially decays to zero, with one small minimum after initially passing through zero. Hence, it is evident that atoms in a liquid do not influence their neighbours as significantly as in the solid phase, which would be expected, owing to the fact that their positions are not fixed and interatomic interactions are slightly weaker.

There is a larger difference between the VACF of the HO and L-J as they behave differently at large displacements. For the L-J, potential energy reduces to zero at large displacements; hence, as atoms part they exert smaller forces and become uncorrelated. Whereas, at large displacements for the HO, force subjected by atoms increase. Therefore, atoms remain within close proximity and have velocities that are perfectly correlated. This is why the VACF files reduce two zero, but the HO does not.

Integration of VACF

Figure 20 displays the velocity autocorrelation function as a function of time for all phases, with 8,000 atoms. As expected, all phases reduce to zero as $\tau$ tends to large values of time. This was expected as velocities at large times should be uncorrelated to initial velocities. Clearly, in the vapour phase, atoms took a long time to become uncorrelated, as seen from the slow decay. Whereas, both liquid and solid phases reduced to zero almost immediately. Hence, the motion of atoms become uncorrelated relatively quickly in these phases; the solid phase becomes uncorrelated before the liquid phase, as expected. Atoms in the vapour phase are disperse, which results in the time taken for interatomic interactions to occur relatively large, despite their large velocities. Solid and liquid phases have atoms that are relatively close together; therefore, interatomic interactions are frequent, which results in the velocities becoming uncorrelated quickly.

Figure 20. VACF as a function of time for solid, liquid and vapour phases, seen in blue, orange and green, respectively.

Figure 21. Running integral approximated by the trapezium rule as a function of time for solid, liquid and vapour phases, seen in blue, orange and green, respectively.

Figure 22. Running integral approximated by the trapezium rule as a function of time for solid, liquid and vapour phases, seen in blue, orange and green, respectively. The range where the solid and liquid phases converge to has been magnified so to allow analysis of their changes with time.

A running integral of the velocity autocorrelation function as a function of time is displayed in figure 21 for 8,000 atoms. It is evident that the vapour phase has significantly larger values than both of the other phases. This suggests that the diffusion coefficient for vapours is significantly larger than that of solids and liquids, which is expected.

Data was also processed for a 1,000,000 atoms, for all three phases, and displayed in figures 23 and 24; the former is the velocity autocorrelation function as a function of time and the latter a running integral of the velocity autocorrelation function as a function of time.

Figure 23. VACF as a function of time for solid, liquid and vapour phases, seen in blue, orange and green, respectively.

Figure 24. Running integral approximated by the trapezium rule as a function of time for solid, liquid and vapour phases, seen in blue, orange and green, respectively.

Figure 25. Running integral approximated by the trapezium rule as a function of time for solid, liquid and vapour phases, seen in blue, orange and green, respectively. The range where the solid and liquid phases converge to has been magnified so to allow analysis of their changes with time.

From figures 21 and 24, the diffusion coefficient was estimated. The table below displays the obtained results for 8,000 and 1,000,000 atoms. There is fairly good agreement between diffusion coefficients for the liquid and vapour phase, so the small number of atoms is sufficient to calculate this parameter for these phases. However, the diffusion coefficient for the solid phase with 8,000 atoms is not of the same order of magnitude as all of the other calculated values.

Phase	D / cm2s-1	D / cm2s-1 (Million)
Solid	$1.33\times 10^{-04}$	$4.55\times 10^{-05}$
Liquid	$9.79\times 10^{-02}$	$9.01\times 10^{-02}$
Vapour	$2.12$	$3.27$

Diffusion coefficient for both calculation methods are in good agreement.

Possible sources of error when calculating the diffusion coefficient were: estimation of the area from the trapezium rule, using a relatively small number of atoms and intrinsic errors to the calculation method. It has already been seen that the number of atoms used is sufficient to obtain accurate results, so this source can not be the largest. Intrinsic errors with the calculation method can only be isolated by comparing the values to experimental or other theories. Hence, this error was neglected. The largest error therefore, is in the trapezium rule. Owing to the large curvature in parts of the velocity autocorrelation functions, estimation of the diffusion coefficient will be smaller than expected.