Add tutorial 1 on modelling a simple polymer chain

tut1_polymer/README

+In this second tutorial, we use LAMMPS to simulate the dynamics of a simple
+polymer (a chain of beads) within a solvent (modelled implicitly by simulating
+the system using Langevin/Brownian dynamics). As in the previous tutorial,
+beads interact with one another via the truncated Lennard-Jones (LJ)
+potential, with the default parameters setting the interactions to be purely
+repulsive.
+
+To link consecutive beads of the polymer together, we add finite extensible
+nonlinear elastic (FENE) springs/bonds between them. In formula, the FENE
+potential is given by
+
+   U(r) = -K_f*R_0^2*ln(1-(r/R_0)^2),
+
+where r is the distance between two consecutive beads, and K_f is the spring
+constant. The FENE bond differs from the harmonic bond in that it can only
+extend to a finite limit (R_0), and this can help reduce the chance of the
+polymer getting entangled.
+
+We also model the polymer to be semi-flexible by adding angle bonds - bonds
+which act between three consecutive beads. Here we use the cosine potential
+
+   U(theta) = K_b*(1-cos(theta)),
+
+where theta is angle formed between the two bonds linking three consecutive
+beads together, and K_b is a parameter that controls the flexibility or the
+persistence length l_p (K_b ~ k_B*T*l_p/sigma), the characteristic length over
+which the polymer loses memory of its direction. This kind of semi-flexible
+polymer models is also known as a Kratky-Porod model.
+
+What to do:
+
+0. Plot the bond and angle potentials
+
+   Using your favourite plotting package, plot both the FENE and cosine
+   potentials. Make sure you are happy with how these potentials do what they
+   are meant to do in the simulations. Explore different parameter values and
+   see how that changes the shape of the potential.
+
+1. Run the simulation with the default parameters
+
+   a. Generate the required input scripts/files
+
+         bash init.sh
+
+      As in the first tutorial, 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 command and each section does! In particular:
+      -  What are some new commands added in the run.lam file?
+      -  What is the initial configuration of the polymer? And how did we
+      	 equilibrate the system?
+      -  What additional info have we provided in the init.in file?
+
+      As always, if in doubt, check out LAMMPS documentation on
+
+         docs.lammps.org/Manual.html
+
+   b. Run the simulation
+
+      This can be done using LAMMPS installed on the school computer
+
+         lmp -in run.lam
+
+      pos*.lammpstrj files will be generated, storing the beads' positions
+      and velocities
+
+   c. Visualise the trajectory of the polymer
+
+      This can be done using VMD as in the first tutorial. Note that because
+      LAMMPS does not output the bead positions in order (i.e., by the
+      bead index) while VMD thinks that is the case, beads are linked together
+      wrongly in the visualisation. We can fix this by typing the following
+      set of commands in the terminal after loading the molecule:
+
+      topo clearbonds
+      set N {nbeads}
+      for {set i 0} {$i < [expr $N-1]} {incr i} {
+      	  set j [expr $i+1]
+	  topo addbond $i $j
+      }
+
+      [You will need to replace the variable nbeads with the actual number of
+       beads in the polymer]
+      
+2. Explore the effect of changing the
+
+   
+
+3. 
+

tut1_polymer/create_polymer.py

+#!/usr/bin/env python3
+# create_polymer.py
+# A simple script to generate a random walk polymer in a box of size (lx,ly,lz)
+# The centre of the box is at the origin and the first bead of the polymer
+# is at point (x0,y0,z0) (or origin by default)
+
+import sys
+import numpy as np
+
+args = sys.argv
+nargs = len(args)
+
+if (nargs != 9 and nargs != 12):
+    print("Usage: create_polymer.py nbeads ntypes sigma lx ly lz seed",
+          "out_file [x0 y0 z0]")
+    sys.exit(1)
+
+nbeads = int(args.pop(1))  # Number of polymer beads
+ntypes = int(args.pop(1))  # Number of bead types
+sigma = float(args.pop(1)) # The bond length between beads (or bead size)
+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
+
+xhalf = lx/2.0
+yhalf = ly/2.0
+zhalf = lz/2.0
+
+xprev = 0.0
+yprev = 0.0
+zprev = 0.0
+
+# A helper function to check if the coordinates lie outside the box
+def out_of_box(x,y,z):
+    global xhalf, yhalf, zhalf
+    return abs(x) > xhalf or abs(y) > yhalf or abs(z) > zhalf
+
+if (nargs == 12):
+    x0 = float(args.pop(1))
+    y0 = float(args.pop(1))
+    z0 = float(args.pop(1))
+    # Check to make sure that the initial point is within the simulation box
+    if (out_of_box(x0,y0,z0)):
+        print("Error: the initial point must be within the simulation box")
+        sys.exit(1)
+    xprev = x0
+    yprev = y0
+    zprev = z0
+
+# Set up the random generator
+rng = np.random.default_rng(seed)
+
+with open(out_file,'w') as writer:
+    # Output the position of the first bead
+    t = rng.integers(1,ntypes+1)
+    writer.write("{:d} {:d} {:f} {:f} {:f}\n".format(1,t,xprev,yprev,zprev))
+    # Generate the positions of the remaining beads
+    for i in range(2,nbeads+1):
+        t = rng.integers(1,ntypes+1)
+        while (True):
+            # Generate a random position for the next bead to be anywhere on a
+            # sphere of diameter sigma that is centred at the previous bead
+            r = rng.random()
+            costheta = 1.0-2.0*r
+            sintheta = np.sqrt(1-costheta*costheta)
+            r = rng.random()
+            phi = 2.0*np.pi*r
+            x = xprev+sigma*sintheta*np.cos(phi)
+            y = yprev+sigma*sintheta*np.sin(phi)
+            z = zprev+sigma*costheta
+
+            # Make sure the new coordinates are not outside the box,
+            # otherwise re-generate the point
+            if (out_of_box(x,y,z)): continue
+            
+            writer.write("{:d} {:d} {:f} {:f} {:f}\n".format(i,t,x,y,z))
+            xprev = x
+            yprev = y
+            zprev = z
+            break
+                

tut1_polymer/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
+nbeads=1000
+ntypes=1
+seed=93845
+
+# Set LAMMPS parameters
+# LJ potential
+# Try changing epsilon and cutoff and see how that affects the dynamics!
+epsilon=1.0
+cutoff=1.122462048309373
+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
+
+# FENE potential
+k_f=30.0
+R_0=1.6
+# These parameter values are used such that the equilibrium bond length
+# is roughly 1.0 sigma
+
+# Angle potential
+# Try changing the persistence
+persist_length=3.0 # Set the persistence length (stiffness) of the polymer
+
+# 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=$(get_norm $epsilon $cutoff $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_1=100
+run_init_time_2=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_1=$($pyx -c "print(int($run_init_time_1/$dt))")
+run_init_time_2=$($pyx -c "print(int($run_init_time_2/$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
+polymer_py="create_polymer.py"
+
+# Generate the lammps input file
+traj_in="traj.in"
+$pyx $polymer_py $nbeads $ntypes $sigma $lx $ly $lz $seed $traj_in
+
+init_file="init.in"
+echo "LAMMPS data file via
+
+${nbeads} atoms
+${ntypes} atom types
+$((${nbeads}-1)) bonds
+1 bond types
+$((${nbeads}-2)) angles
+1 angle 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 # angle"
+print ""
+}{
+print NR,1,$2,$3,$4,$5,0,0,0
+} END {
+print ""
+print "Bonds"
+print ""
+for (i = 1; i <= NR-1; i++) {
+print i,1,i,i+1
+}
+print ""
+print "Angles"
+print ""
+for (i = 1; i <= NR-2; i++) {
+print i,1,i,i+1,i+2
+}}' $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 angle # State the system includes bonds and 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
+
+# Commands setting how MPI processes communicate with each other, should one
+# choose to run the simulation in parallel using MPI. They are used here to
+# avoid a warning saying that the communication cutoff length is smaller than
+# the bond length
+comm_style tiled
+comm_modify mode single cutoff 4.0 vel yes
+
+# 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)
+compute gyr all gyration # Compute the radius of gyration of the polymer
+thermo_style custom step temp epair emol vol c_gyr
+
+# 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 harmonic bonds to link between beads 
+# U_harm(r) = k*(r-r_0)
+bond_style harmonic
+# Key params: bond_type k r_0
+bond_coeff 1 100.0 1.1
+
+# Use cosine angle bonds to model a semi-flexible polymer
+# U(theta) = K_b*(1-cos(theta))
+angle_style cosine
+# Key params: angle_type K_b
+angle_coeff 1 10.0 # Make the fibre very stiff first to help remove overlap
+
+# Use the soft potential for equilibration to push apart atoms close together
+pair_style soft ${cutoff_0}
+pair_coeff * * 100.0 ${cutoff_0}
+
+##################################################
+
+# Set integrator/dynamics
+
+# Use the NVE ensemble (particle number, volume, energy conserved)
+fix 1 all nve
+# Use the langevin thermostat for performing 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_1}
+
+##################################################
+
+# Equilibrate with FENE bonds
+
+# Switch from harmonic to FENE bonds now that the polymer has relaxed a bit
+bond_style fene
+# Key params: k_f R_0 epsilon sigma
+bond_coeff 1 ${k_f} ${R_0} 1.0 ${sigma}
+
+# The following command is needed so as to switch off the pairwise (1-2)
+# interactions between consecutive beads that are set from using the
+# pair_style command. This is because the implementation of FENE bonds in 
+# LAMMPS already includes a WCA repulsive part, and including the 1-2 
+# interactions between these bonded beads will lead to double-counting.
+# See LAMMPS documentation for more details
+special_bonds fene
+
+# Change the polymer to the desired stiffness
+angle_coeff 1 ${persist_length}
+
+# Switch from soft to the purely repulsive LJ (or WCA) potential
+pair_style lj/cut ${cutoff_0}
+pair_coeff * * 1.0 ${sigma} ${cutoff_0}
+
+run ${run_init_time_2}
+
+##################################################
+
+# Main simulation
+
+# Change the LJ interactions from purely repulsive to the desired strength
+pair_style lj/cut ${cutoff_0}
+pair_coeff * * ${epsilon} ${sigma} ${cutoff}
+
+##################################################
+
+# 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