Add tutorial 0 on simulating a simple LJ fluid

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.
+
+What to do:
+
+1. Run the simulation to get a feel of how to use LAMMPS
+
+   a. Generate the required input scripts/files
+
+      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)
+
+      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
+      linux computers. Alternatively, you can also compile LAMMPS on your
+      own computer. Visit docs.lammps.org/Install.html for more info on how to
+      install LAMMPS. 
+
+      lmp -in *.lam
+
+      pos*.lammpstrj files will be generated storing the atoms' positions
+      and velocities. See the 'dump' command in the *.lam script.
+
+   c. Visualise the atoms' trajectories
+
+      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
+      
+      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]
+	
+      - 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:
+
+   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

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
+
+import sys
+import numpy as np
+
+args = sys.argv
+
+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)
+
+rng = np.random.default_rng(seed)
+
+xhalf = lx/2.0
+yhalf = ly/2.0
+zhalf = lz/2.0
+
+with open(out_file,'w') as writer:
+    for i in range(1,natoms+1):
+        t = rng.integers(1,ntypes+1)
+        x = rng.random()*lx-xhalf
+        y = rng.random()*ly-yhalf
+        z = rng.random()*lz-zhalf
+        writer.write("{:d} {:d} {:f} {:f} {:f}\n".format(i,t,x,y,z))

tut0_lj_fluid/init.sh

+#!/bin/bash
+# init.sh
+# A script to generate the lammps script and init config file
+
+# Make sure we are using python3
+pyx=python3
+
+# Set the box size, number of atoms and random seed
+lx=50
+ly=50
+lz=50
+natoms=1000
+ntypes=2
+seed=3432
+
+# Set LAMMPS parameters
+# LJ potential
+# Try changing epsilon and cutoff and see how that affects the dynamics!
+epsilon_11=1.0
+epsilon_12=1.0
+epsilon_22=1.0
+cutoff_11=1.8
+cutoff_12=1.8 #1.122462048309373
+cutoff_22=1.8
+cutoff_0=1.122462048309373 # 2^(1/6) - the cutoff distance for a purely
+                           # repulsive (WCA) potential
+
+sigma=1.0 # Make all atoms have the same size
+
+# Rescale epsilon such that the minimum of the truncated LJ potential
+# actually reaches -epsilon
+function get_norm() {
+    local e=$1
+    local rc=$2
+    local s=$3
+    echo $($pyx -c "
+norm = 1.0+4.0*(($s/$rc)**12.0-($s/$rc)**6.0)
+print($e/norm if norm > 0.0 else $e)")
+}
+epsilon_11=$(get_norm $epsilon_11 $cutoff_11 $sigma)
+epsilon_12=$(get_norm $epsilon_12 $cutoff_12 $sigma)
+epsilon_22=$(get_norm $epsilon_22 $cutoff_22 $sigma)
+
+# Dump frequency and length of simulation
+dt=0.01 # Set timestep size in LJ simulation time unit
+dump_printfreq=100 # In LJ simulation time unit
+thermo_printfreq=100
+run_init_time=100
+run_time=10000
+
+# Convert from simulation time unit to timesteps
+dump_printfreq=$($pyx -c "print(int($dump_printfreq/$dt))")
+thermo_printfreq=$($pyx -c "print(int($thermo_printfreq/$dt))")
+run_init_time=$($pyx -c "print(int($run_init_time/$dt))")
+run_time=$($pyx -c "print(int($run_time/$dt))")
+
+# Set the box boundaries
+xlo=$($pyx -c "print(-int(${lx}/2.0))")
+xhi=$($pyx -c "print(int(${lx}/2.0))")
+ylo=$($pyx -c "print(-int(${ly}/2.0))")
+yhi=$($pyx -c "print(int(${ly}/2.0))")
+zlo=$($pyx -c "print(-int(${lz}/2.0))")
+zhi=$($pyx -c "print(int(${lz}/2.0))")
+
+# Output files
+pos_equil_file="pos_equil.lammpstrj"
+pos_file="pos.lammpstrj"
+
+# Required programs and files
+atom_py="create_atoms.py"
+
+# Generate the lammps input file
+traj_in="traj.in"
+$pyx $atom_py $natoms $ntypes $lx $ly $lz $seed $traj_in
+
+init_file="init.in"
+echo "LAMMPS data file via
+
+${natoms} atoms
+${ntypes} atom types
+
+${xlo} ${xhi} xlo xhi
+${ylo} ${yhi} ylo yhi
+${zlo} ${zhi} zlo zhi
+" > $init_file
+
+awk -v nt=$ntypes 'BEGIN {
+print "Masses"
+print ""
+for (i = 1; i <= nt; i++) {
+print i,1
+}
+print ""
+print "Atoms"
+print ""
+}{
+print NR,$2,$3,$4,$5,0,0,0
+}' $traj_in >> $init_file
+rm $traj_in
+
+# Generate the script to run lammps
+lam_file="run.lam"
+echo "
+##################################################
+
+# Simulation basic setup
+
+units lj # Set simulation unit - for LJ, distance = sigma, energy = k_bT
+atom_style atomic # State the system is only made of atoms without bonds/angles
+boundary p p p # Set periodic boundary conditions
+
+neighbor 1.9 bin # Set how lammps creates neighbour lists
+neigh_modify every 1 delay 1 check yes # How often neighbour lists are created
+
+# Read the file with atoms' init positions
+read_data $(basename $init_file)
+
+##################################################
+
+# Dumps
+
+# Set the params for dumping thermodynamic properties of the system to
+# screen or log file
+thermo ${thermo_printfreq} # Set the output frequency
+# Set the specific observables to dump (see lammps manual for more details)
+thermo_style custom step temp epair emol vol 
+
+# Set up a dump to output atoms' positions
+# xs,ys,zs - positions normalised by the box size (so between 0.0 and 1.0)
+# ix,iy,iz - number of times the atoms crossed each periodic boundary
+dump 1 all custom ${dump_printfreq} $(basename $pos_equil_file) &
+id type xs ys zs ix iy iz
+
+##################################################
+
+# Potentials
+
+# Use the repulsive soft potential for equilibration to push apart atoms that 
+# are close together and would cause numerical divergence had one applied the
+# LJ potential directly
+# U(r) = A*(1+cos(pi*r/r_c))
+# A = strength of repulsion; r_c = cutoff
+pair_style soft ${cutoff_0}
+
+# Set the params for the pairwise interactions between individual atom types
+# Key params: atom_type_1 atom_type_2 A r_c
+pair_coeff * * 100.0 ${cutoff_0} # This sets all pairwise interactions
+
+##################################################
+
+# Set integrator/dynamics
+
+# Use the NVE ensemble (particle number, volume, energy conserved)
+fix 1 all nve
+# Use langevin thermostat for Brownian dynamics
+fix 2 all langevin 1.0 1.0 1.0 ${seed} 
+
+##################################################
+
+# Initial equilibration
+
+# Set the timestep and run the simulation
+timestep ${dt}
+run ${run_init_time}
+
+##################################################
+
+# Main simulation
+
+# Use a (truncated) LJ potential for the interactions between atoms
+# U_LJ(r) = 4.0*epsilon*((sigma/r)^12-(sigma/r)^6)
+# U_LJ/cut(r) = U_LJ(r)-U_LJ(r_c) where r_c is the location of the cutoff
+# If r_c = 2^(1/6) ~ 1.12246, this becomes as purely repulsive (WCA) potential
+# If r_c > 2^(1/6), then there is a region where the potential is attractive
+pair_style lj/cut ${cutoff_0}
+
+# Set the params for the pairwise interactions between individual atom types
+# Key params: atom_type_1 atom_type_2 epsilon sigma cutoff
+# Try changing epsilon and cutoff and see how that affects the dynamics!
+pair_coeff 1 1 ${epsilon_11} ${sigma} ${cutoff_11}
+pair_coeff 1 2 ${epsilon_12} ${sigma} ${cutoff_12}
+pair_coeff 2 2 ${epsilon_22} ${sigma} ${cutoff_22}
+
+##################################################
+
+# Dumps (same settings as above)
+
+undump 1 # Unset previous dump settings
+dump 1 all custom ${dump_printfreq} $(basename $pos_file) &
+id type xs ys zs ix iy iz
+
+##################################################
+
+# Reset time and run the main simulation
+reset_timestep 0
+run ${run_time}
+" > $lam_file