Abstract: Molecular dynamic simulations with 1000 atoms were performed. For the simulations, the
program espresso was used and VMD was used for the visual representation of the atoms. In a first
simulation the movement of the atoms at constant volume could be observed and a melting temperature
of 𝑇𝑚∗ = 1.2 and a boiling temperature of 𝑇𝑚∗ = 1.7 could be estimated. For a second simulation a melting
temperature of 𝑇𝑚∗ = 1.30 and a boiling temperature of 𝑇𝑏∗ = 2.44 could be determined at constant
pressure and 𝜀0 = 2.0. At 𝜀0 = 3.0 a melting temperature of 𝑇𝑚∗ = 2.23 and a boiling temperature of
𝑇𝑏∗ = 3.52 could be determined. Using pair correlation functions, the number of adjacent atoms in the
solid, liquid and gas phase could be determined. These were 𝑛 = 12, 𝑛 = 9 and 𝑛 = 2 respectively.
Thanks to computer simulations, theory and experiments in the laboratory can be complemented. It
is possible to understand the processes in nature and technology more easily. For example,
experiments that are too toxic or dangerous can be simulated by computers. There are several ways
to simulate physical systems. In Molecular Dynamics (MD), the motions of atoms are described using
an equation of motion and one or more potentials acting between those particles. According to
molecular dynamics the time average of the observables is determined. In contrast, in Monte Carlo
simulation, a finite number of random experiments provides the mean value of the observables. In
MD, both momentum and coordinates are attributed to atoms. The trajectories of the atoms can be
determined by using the equation of motion of the Lagrange formalism (1):
d 𝜕𝜕𝓛 𝜕𝓛
( ) − ( ) = 0. (1)
d𝑡 𝜕𝐫̇𝑖 𝜕𝐫𝑖
Here 𝜕𝓛 stands for the Lagrange function, 𝐫𝑖 for the location of the particle 𝑖 and 𝐫̇𝑖 for the time
derivative of the location 𝐫𝑖 . The Lagrange function consists of a potential 𝑉(𝐫𝑖 ) and kinetic term 𝑇(𝐫̇ ):
𝓛(𝐫̇𝒊 , 𝐫𝒊 ) = 𝑇(𝐫̇𝑖 ) − 𝑉(𝐫𝑖 ) . (2)
1 1
𝑇(𝐫̇𝑖 ) = 𝑚𝑖 𝐯𝑖2 = 𝑚𝑖 𝐫̇𝑖2 . (3)
2 2
𝛻𝑉(𝐫𝑖 ) = −𝐅𝑖 . (4)
In equation (3) 𝑚𝑖 stands for the mass and 𝑣𝑖 stands for the velocity of particle 𝑖. In equation (4) 𝐹𝑖
stands for the force acting on a particle 𝑖. By inserting equation (2) in (1) the Newtonian equation of
motion is obtained:
𝑚𝑖 𝐫̈𝑖 = 𝐅𝑖 = −𝛻𝑉𝑖 . (5)
To solve this equation the Velocity-Verlet algorithm can be used. It is used to calculate the positions
1
𝐫𝑖 (𝑡 + 𝛿𝑡), the velocities 𝐯𝑖 (𝑡 + 𝛿𝑡), the accelerations 𝐚𝑖 (𝑡 + 𝛿𝑡) using the Lennard-Jones-Potential
2
and the new velocities 𝐯𝑖 (𝑡 + 𝛿𝑡). These calculations are performed until the simulation time expired.
2.2 Lennard-Jones Potential
The potential energy between molecules or atoms can be described by the Lennard-Jones
Potential 𝑉(𝑟𝑖𝑗 ).
12 6
𝜎0 𝜎0
𝑉(𝑟𝑖𝑗 ) = 4𝜀0 [( ) −( ) ] (7)
𝑟𝑖𝑗 𝑟𝑖𝑗
Here 𝜎0 is the effective diameter, 𝑟𝑖𝑗 the distance between particle 𝑖 and 𝑗 and 𝜀0 is the depth of the
potential well. The first term in brackets indicates the repulsion force and the second the attraction
force. Figure 1 shows a schematic curve for the potential between two molecules.
2
,The potential energy decreases the closer the molecules get to each other. In the ground state distance
𝑟𝑒 , the potential energy has its minimum. If the molecules move closer and closer together, the
dominant force is the repulsion force and the potential energy increases drastically [2].
Figure 1: The schematic potential curve of a Lennard-Jones potential 𝑉(𝑟𝑖𝑗 ) of two particles 𝑖 and 𝑗 depending on particle
distance 𝑟𝑖𝑗 [1].
2.3 Reduced Properties
The introduction of reduced properties allows the description of the thermodynamic state of the
system. For a system containing only identical particles, the fundamental units such as the mass 𝑚,
the potential well depth 𝜀0 as a unit of energy and the minimum distance between two particles 𝜎0 as
a unit of length are set to 1. Thanks to this simplification all other observables can be described by
them. These thermodynamic observables are recorded in the following table [1]:
Table 1: Reduced dimensions based on the fundamental units of mass 𝑚, of energy 𝜀0 and of length 𝜎0 [1].
Density 𝝆∗ = 𝝆𝝈𝟑𝟎
Temperature 𝑇 ∗ = 𝑘𝐵 𝑇/𝜖0
Energy 𝐸 ∗ = 𝐸/𝜖0
Pressure 𝑝∗ = 𝑝𝜎03 /𝜖0
Time 𝑡 ∗ = (𝜖0 /𝑚𝜖0 )1/2𝑡
Force 𝐟 ∗ = 𝐟𝜎0 /𝜖0
2.4 Ensembles and Ergoden Hypothesis
In statistical thermodynamics, macroscopic systems, such as temperature, are studied from a
molecular point of view, such as the kinetic energy of individual atoms. A macroscopic system can be
described by one or more microstates. These are defined by the stationary quantum states 𝑖 with
associated wave function 𝜓𝑖 and energy eigenvalue 𝐸𝑖 . A distinction is made between two different
many-body systems. In the so-called microcanonical ensemble, the number of particles 𝑁, the volume
𝑉 and the energy 𝐸 are kept constant. In the canonical ensemble, however, the number of particles 𝑁,
volume 𝑉 and the temperature 𝑇 are kept constant. In the microcanonical ensemble, microstates of
equal energies are considered, and an average value over time is obtained. According to the quasi-
ergodic hypothesis, however, it does not matter whether the microstates are observed as the average
over time or as the average over a certain number of copies of a system [1].
3
,2.5 Pair Correlation Function
The pair correlation function 𝑔(𝑟) is used to determine the short- and long-range order of a fluid.
𝜌(𝑟)
𝑔(𝑟) = (8)
𝜌̅
The average particle density of a fluid is determined by the following equation:
𝑁
𝜌̅ = (9)
𝑉
The crystalline phase has both short- and long-range order, whereas the liquid phase loses its long-
range order. In the gaseous phase there is neither short- nor long-range order. This is because the
greater disorder in the gas phase means that the information about the short- and long-range order is
lost.
4
,2 Experiment
2.1 Task
The behavior of 1000 atoms during heating is to be simulated with the help of MD simulation. In the
first simulation the volume is kept constant. In two others the pressure is kept constant for
𝜀0 = 2.0 and 𝜀0 = 3.0. The pair correlation functions 𝑔(𝑟) at different temperatures should be
determined and discussed. In addition, for each phase the local density function is to be integrated
over the volume of the first coordination shell and the number of adjacent atoms is to be calculated.
2.2 Implementation
For the simulations, the Velocity-Verlet algorithm and the Lennard-Jones potential for the interatomic
interactions were used. The parameters used in the respective simulations are recorded in the
following table:
In the first simulation, the phase transition temperatures were determined by observation of the
atoms. You can see from the beginning of the simulation that the atoms are still in the solid phase.
They oscillate but are still fixed in their lattice positions. A phase transition to the liquid phase can be
recognized at 𝑇 = 1.2 by the fact that the atoms now move faster and are no longer fixed at their
positions. At 𝑇 = 1.7 the atoms are now moving even faster, and the atoms are thus now in the gas
phase.
3.2 Simulation at constant pressure
To determine the melting and boiling temperatures, the volume was plotted against the temperature.
Figure 2 shows the graph for the potential depth at 𝜀0 = 2.0 and Figure 3 shows the graph at 𝜀0 = 3.0.
The solid phase of the systems extends to the first leap point of the graph. This data was extrapolated
by linear regression (line of best fit 𝑦1 ). After the first leap, the system changes to the liquid phase. A
linear regression was also performed here (line of best fit 𝑦2 ). The volume hardly increases at all, since
only the short-range order is lost, but not the long-range order. After the second leap in the graph the
change to the gaseous phase takes place. Here, linear extrapolation was also performed (line of best
fit 𝑦3 ). It can also be seen that now the short-range order has also been lost. This is noticeable by the
strong increase in volume. The reduced melting temperature 𝑇𝑚∗ is determined by the point of
intersection of 𝑦1 and 𝑦2 :
704.17 𝑇𝑚∗ + 159.93 = 119.46 𝑇𝑚∗ + 920.27 (10)
920.27 − 159.93
𝑇𝑚∗ = = 1.3
704.17 − 119.46
The reduced boiling temperature 𝑇𝑚∗ is determined by the point of intersection of 𝑦2 and 𝑦3 . The
calculations are carried out in the same way for 𝜀0 = 3.0 and summarized in Table 3.
Table 3: The reduced melting 𝑇∗𝑚 and boiling temperature 𝑇∗𝑏 of the different pontetial depths.
𝜺𝟎 𝑻∗𝒎 𝑻∗𝒃
2.0 1.30 2.44
3.0 2.23 3.52
A comparison of the data from Table 3 shows that the determined temperatures are higher for 𝜀0 = 3.0
than for 𝜀0 = 2.0. This result makes sense, because at a greater potential well the interaction between
the particles is greater. Therefore, the system is more stable, and more energy is needed to enable the
phase transitions.
6
The benefits of buying summaries with Stuvia:
Guaranteed quality through customer reviews
Stuvia customers have reviewed more than 700,000 summaries. This how you know that you are buying the best documents.
Quick and easy check-out
You can quickly pay through credit card or Stuvia-credit for the summaries. There is no membership needed.
Focus on what matters
Your fellow students write the study notes themselves, which is why the documents are always reliable and up-to-date. This ensures you quickly get to the core!
Frequently asked questions
What do I get when I buy this document?
You get a PDF, available immediately after your purchase. The purchased document is accessible anytime, anywhere and indefinitely through your profile.
Satisfaction guarantee: how does it work?
Our satisfaction guarantee ensures that you always find a study document that suits you well. You fill out a form, and our customer service team takes care of the rest.
Who am I buying these notes from?
Stuvia is a marketplace, so you are not buying this document from us, but from seller DevChem. Stuvia facilitates payment to the seller.
Will I be stuck with a subscription?
No, you only buy these notes for $11.27. You're not tied to anything after your purchase.