Table of Contents

Lennard-Jones liquids

In this exercise, you will simulate a fluid of monoatomic particles that interact with a Lennard-Jones potential. The method to be used is molecular dynamics (MD) with periodic boundary conditions using CP2K. The aim is to explore the method, calculate the
radial distribution function $g(r)$

You are expected to hand in the respective plots by email, ONLY PDF format.

Run your first simulation using CP2K

When you check CP2K's features and the outline of the lecture you will notice that there are many levels of theory, methods, and possibilities to combine them. This results in a large number of possible options and coefficients to setup, control and tune a specific simulation.
Together with the parameters for the numerical solvers, this means that an average CP2K configuration file will contain quite a number of options (even though for many others the default value will be applied) and not all of them will be discussed in the lecture or the exercises.

The CP2K Manual is the complete reference for all configuration options. Where appropriate you will find a reference to the respective paper when looking up a specific keyword/option.

To get you started, we will do a simple exercise using Molecular Mechanics (that is: a classical approach). The point is to get familiar with the options, organizing and editing the input file and analyze the output.

Background

You are expected to carry out an MD simulation of Lennard-Jones fluid containing mono-atomic
particles. The potential energy expression used is

where is $\epsilon$ the well depth, σ is related to the minimum of Lennard-Jones potential, rij is the
distance between atoms i and j.

Periodic boundary conditions should be used in this simulation. The atom near the ”edge” of the simulation box interacts with atoms contained in the neighboring image of the box. In computer simulations, one of these is the original simulation box, and others are copies called images. During the simulation, only the properties of the original simulation box need to be recorded and propagated. The minimum-image convention is a common form of PBC particle bookkeeping in which each individual particle in the simulation interacts with the closest image of the remaining particles in the system. To prevent artifacts, it requires a cut-off value of rij for the L-J potential. For realistic results, the cut-off should be less than the half of the simulation box size and over σ Radial distribution
function, (or pair correlation function) $g(r)$ in a system of particles (atoms, molecules, colloids,
etc.), describes how density varies as a function of distance from a reference particle.

Part I: Set up MD simulation

In this section, a commented CP2K input example for an MD calculation is provided.
Comments are added, they start with a hash mark '#'.

1. Step

Load the CP2K module as shown before, create a directory ex1 and change to it:

mkdir ex1
cd ex1

Save the following input to a file named argon.inp (for example using $ vim argon.inp):

Instructions starting with an ampersand & start a section and must be terminated with an &END SECTION-NAME. Other instructions are called keywords. The indentation is ignored but recommended for readability.

2. Step

Run a CP2K calculation with the following command:

$ cp2k.sopt -i argon.inp -o argon.out

Alternatively you can run it using

$ cp2k.sopt -i argon.inp | tee argon.out

which will write the output simultaneously to a file argon.out and show it in the terminal.

TASK

Run the calculation and visualize the trajectories using VMD

Run calculations with different timesteps (0.5 2 5fs), different temperatures(84, 300, 400K), different densities($\rho$=0.25,0.5,1), and a different number of timesteps (100,1000,5000), and inspect geometry in each case, plot the total energies, temperature, potential energies. Try to comment and explain what you observe.

Part II: Force Field Parameter

You need to modify the L-J force field parameter. Run calculations with different $\sigma$ and $\epsilon$ and inspect the insight of L-J potential (Energies, velocities, temperatures etc ). Then use the force field parameter provided to plot the L-J potential curve.

… and using > energy_${d}A.out we redirect the output of the sed command to new files energy_2.0A.out, energy_2.1A.out, etc.

Then we run cp2k.sopt as shown before on those new input files and write the output to new output files as well

Using awk we extract the energy from the output file

TASK

Plot the Lennard-Jones potential against Ar-Ar distance $r$ (2-4 Å) using different $\epsilon$ and $\sigma$.

Repeat the L-J MD calculation with different $\epsilon$ and $\sigma$, and compare the potential energy, temperature evolution, and explain the relation between calculated properties and force field parameters.

Part III: Radial distribution functions

Use VMD or write your own program (Fortran, C, C++, Python etc.) to calculate radial distribution $g(r)$. Plot $g(r)$, and against various the temperatures to examine the effects.
VMD comes with an extension for exactly this purpose: In the VMD Main window open “Extensions → Analysis” click on “Radial Pair Distribution function $g(r)$. In the appearing window use “Utilities → Set unit cell dimensions” to let VMD know the simulation box you used. After that use Selection 1 and 2 to define the atomic types that you want to calculate the rdf for, for example “element Ar”. In the plot window, use “File”, you can save the plot data.

TASK

Plot $g(r)$ at 84, 300 and 400 K into the same graph.

What are the differences in the height of the first peak, and why does temperature contribute to the differences?

Compared to experimental data exp_gr.dat taken at 84 K, what does this say about the structure of the liquid and is this expected? .