Add more contents and comments to tutorial 0

tut0_lj_fluid/README

-This is a simple example of using LAMMPS to perform molecular dynamics
-simulations of a Lennard-Jones (LJ) fluid. Atoms are simulated in a periodic
-box, and they interact with each other via the (truncated) LJ potential.
+In this first tutorial, we use LAMMPS to perform molecular dynamics (MD)
+simulations of a Lennard-Jones (LJ) fluid (colloidal particles in a sovlent).
+The particles (atoms) are simulated in a periodic box, and they interact with
+each other via the truncated LJ potential, which is given by
+
+   U_LJ/cut(r) = U_LJ(r)-U_LJ(r_c), if r < r_c
+
+and 0 otherwise. Here, r is the distance between two particles, r_c is the
+cutoff distance at which point the potential is shifted and truncated, and
+U_LJ is the usual Lennard-Jones potential
+
+   U_LJ(r) = 4*epsilon*((sigma/r)^12-(sigma/r)^6),
+
+with epsilon being the interaction strength and sigma the distance where U = 0
+(the bead size).
+
+Rather than simulating all the solvent particles (which is computationally
+very expensive!), we model their effect on the colloidal particles implicitly
+by including more terms to the equations of motion, in addition to conservative
+forces (F_con), such that
+
+   m*dv(t)/dt = F_con + F_damp + F_rand
+
+with
+
+   F_con = -nabla(U)
+   F_damp = -gamma*v(t)
+   F_rand = sqrt(2.0*gamma*k_B*T)*eta(t).
+
+Here, F_damp is a damping force on the particle (with gamma being the damping
+const and v the particle's velocity), and F_rand is a random force on the
+particle due to collisions with the solvent particles. eta(t) is a random
+'noise' vector that obeys the following statistical averages
+
+   <eta(t)> = 0
+   <eta_i(t)eta_j(t')> = delta(t-t')delta_ij
+
+where the first delta in the second equation is a dirac delta and the second
+one is a Kronecker delta, and i,j are indices running over Cartesian
+components. This type of equations is known as Langevin equations, and you
+will learn more about them in the Year 4 Statstical Physics course. Often MD
+simulations using this scheme are referred as 'Brownian/Langevin dynamics'
+simulations.
+
+If you are unsure at any point of these tutorial exercises, please consult
+the online LAMMPS documentation:
+
+   docs.lammps.org/Manual.html
+
+==============================================================================
 
 What to do:
 
+0. Familarise yourself with the truncated LJ potential
+
+   Using your favourite plotting package, plot the truncated LJ potential for
+   different parameter values (epsilon and r_c) and have a think about the
+   following questions:
+   -  At what critical value for r_c does the potential become purely
+      repulsive? Note that setting r_c to this value (with epsilon = 1.0)
+      gives a potential also known as the Weeks-Chandler-Anderson (WCA)
+      potential. [Hint: have a look at the provided bash script init.sh if you
+      are unsure.]
+   -  How does shifting the LJ potential change the value of its global
+      minimum? How does this depend on r_c? Is there a way we can ensure that
+      the minimum remains the same as that for the non-truncated LJ potential
+      (i.e., that it reaches -epsilon)?
+	 
 1. Run the simulation to get a feel of how to use LAMMPS
 
    a. Generate the required input scripts/files
 
-      bash init.sh
+         bash init.sh
 
       This generates two files that are necessary for running LAMMPS
-      - run.lam file - a 'driver' script that is read by LAMMPS
-      - init.in file - a trajectory file containing the initial configuration
-                     of the system (i.e., the atoms' positions, bonds, etc)
+      -  run.lam file - a 'driver' script that is read by LAMMPS
+      -  init.in file - a trajectory file containing the initial configuration
+                        of the system (i.e., the atoms' positions, bonds, etc)
 
       Have a look at these files individually and make sure you understand
       what each function and each section does!
 
-      If in doubt, check out LAMMPS documentation on
-
-      docs.lammps.org/Manual.html
-
    b. Run the simulation
 
       You can use the pre-compiled version of LAMMPS available on the school
@@ -29,7 +87,7 @@ What to do:
       own computer. Visit docs.lammps.org/Install.html for more info on how to
       install LAMMPS. 
 
-      lmp -in *.lam
+         lmp -in *.lam
 
       pos*.lammpstrj files will be generated storing the atoms' positions
       and velocities. See the 'dump' command in the *.lam script.
@@ -39,40 +97,64 @@ What to do:
       This can be done using VMD on the school computers. You can also
       download VMD using this link:
 
-      https://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD
+      	 www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD
       
       You can visualise the trajectory file as follows:
 
-      - Type 'vmd' in the terminal
-      - In the VMD Main window, click on 'File' -> 'New Molecule...'. In the
-      	pop-up window, click on 'Browse...' and select the trajectory file
-	(pos.lammpstrj). Click 'Load' to load the data
-      - The display will load the unwrapped positions of the atoms. To 'rewrap'
-      	their positions, we need to type the following in the terminal:
-
-	pbc set {lx ly lz} -all
-	pbc wrap -all
-
-	[You will need to replace lx ly lz with the actual box size]
+      -  Type 'vmd' in the terminal
+      -  In the VMD Main window, click on 'File' -> 'New Molecule...'. In the
+      	 pop-up window, click on 'Browse...' and select the trajectory file
+	 (pos.lammpstrj). Click 'Load' to load the data
+      -  The display will load the unwrapped positions of the atoms. To
+      	 'rewrap' their positions, we need to type the following in the
+	 terminal:
+	 
+	 pbc set {lx ly lz} -all
+	 pbc wrap -center origin -all
+
+	 [You will need to replace lx ly lz with the actual box size]
 	
-      - To get a nicer representation, you can make the following changes. In
-      	the VMD Main window, click on 'Graphics' -> 'Representations...'.
-      	In the pop-up window, select 'VDW' as the Drawing Method and use
-	'AOChalky' for the Material. Feel free to play around with other
-	display settings!
+      -  To get a nicer representation, you can make the following changes. In
+      	 the VMD Main window, click on 'Graphics' -> 'Representations...'. In
+	 the pop-up window, select 'VDW' as the Drawing Method and use
+	 'AOChalky' for the Material. Feel free to play around with other
+	 display settings!
 
 2. Play around with the simulation parameters
 
    Explore how changing the number (density) of atoms and the interaction
-   strength between different types of atoms change the overall dynamics
-   of the system. This can be done by editing the parameter values in the
-   init.sh script, in particular:
+   strength between different types of atoms changes the overall dynamics of
+   the system. This can be done by editing the parameter values in the init.sh
+   script, in particular
 
    natoms     = number of atoms in total
    epsilon_ab = the interaction strength between type a and type b atoms
    cutoff_ab  = the cutoff distance of the potential for the interaction
                 between type a and type b atoms
 
-   See if you can find a range of parameter values where:
-   - the two types of atoms are well mixed
-   - the two types of atoms form their own clusters
\ No newline at end of file
+   See if you can find a range of parameter values where
+   -  the two types of atoms are well mixed
+   -  the two types of atoms form their own clusters
+
+   For a more interesting challenge, try modifying init.sh to simulate a
+   system with 3 (or more) types of atoms!
+
+3. Understand the units used in the simulations
+
+   Simulations are often conducted in 'reduced' units to maintain generality
+   of the results and to limit errors introduced from floating-point
+   arithmetic caused by doing computation on variables spanning many different
+   orders of magnitude. This is achieved by rescaling variables/observables by
+   a set of fundamental quantities characteristic to the system. In the case
+   of an LJ fluid, three key properties are the energy, length, and mass of
+   the system, and a natural choice is to rescale these by the thermal energy
+   (k_B*T), the size of a colloidal particle (sigma), and its mass (m). In
+   other words, we have k_B*T = 1.0, sigma = 1.0 and m = 1.0 in
+   reduced/simulation units (also known as LJ units for this particular
+   rescaling in LAMMPS).
+   -  Based on these fundamental units, can you construct a natural time
+      scale to describe the system? [Hint: take a look at the 'units' command
+      in LAMMPS documentation.]
+   -  A more challenging task - looking at the dimensions of the variables
+      involved in the equations of motion above, can you identify other
+      relevant time scales for describing the system?

tut0_lj_fluid/create_atoms.py

 #!/usr/bin/env python3
 # create_atoms.py
-# A simple script to generate atoms in a box of size (lx,ly,lz)
-# The centre of the box is at the origin
+# A simple script to generate atoms in a box of size (lx,ly,lz). The centre of
+# the box is at the origin
 
 import sys
 import numpy as np
@@ -12,20 +12,22 @@ if (len(args) != 8):
     print("Usage: create_atoms.py natoms ntypes lx ly lz seed out_file")
     sys.exit(1)
 
-natoms = int(args.pop(1))
-ntypes = int(args.pop(1))
-lx = float(args.pop(1))
-ly = float(args.pop(1))
-lz = float(args.pop(1))
-seed = int(args.pop(1))
-out_file = args.pop(1)
+natoms = int(args.pop(1)) # Number of atoms
+ntypes = int(args.pop(1)) # Number of atom types
+lx = float(args.pop(1))   # Box dimension in the x direction
+ly = float(args.pop(1))   # Box dimension in the y direction
+lz = float(args.pop(1))   # Box dimension in the z direction
+seed = int(args.pop(1))   # Seed for the random generator
+out_file = args.pop(1)    # Name of the output file
 
+# Set up the random generator
 rng = np.random.default_rng(seed)
 
 xhalf = lx/2.0
 yhalf = ly/2.0
 zhalf = lz/2.0
 
+# Generate and output the atoms' positions
 with open(out_file,'w') as writer:
     for i in range(1,natoms+1):
         t = rng.integers(1,ntypes+1)