LAMMPS is open-source Simulation Package developed at Sandia National Laboratories for material designing.

How to install LAMMPS using a pre-built executable

Use the following commands in a Linux terminal:

1. $ sudo add-apt-repository ppa:gladky-anton/lammps
2. $ sudo add-apt-repository ppa:openkim/latest
3. $ sudo apt-get update
4. $ sudo apt-get install lammps-stable

How to build LAMMPS with extra-packages using cmake

1. Download the recent stable release from the website as a tarball. The latest release is LAMMPS Stable Release-23rd June.
2. Install Cmake ($ sudo apt install cmake) if you havn't already.
3. Extract the downloaded tar.gz file. ($ tar -xf lammps-stable.tar.gz)
4. create a directory called build. ($ mkdir build)
5. Go inside the build directory. ($ cd build)
6. Type $ cmake -C ../lammps-23Jun2022/cmake/presets/most.cmake -C ../lammps-23Jun2022/cmake
7. Optionally, install additional packages (Example: LATBOLTZ package) using ($ cmake -D PKG_LATBOLTZ=ON .)
8. Finally, type "$ make install" which will generate an executable file named lmp.

Repulsive Lennard Jones Particle Simulation

Create an input file: You need to create a data file that contains the specifications of your particles. The data file should include the following information:

1. The number of particles
2. The mass of each particle
3. The initial positions and velocities of each particle

# Specify simulation box size and shape
units           lj
dimension       3
boundary        p p p
atom_style      atomic

# Specify the number of particles and their initial positions and velocities
read_data       particles.data

# Specify the potential energy function
pair_style      lj/cut 2.5
pair_coeff      * * 1.0 1.0 2.5

# Specify the timestep and total simulation time
timestep        0.005
run             1000

 
# Lennard-Jones particle data

# Specify the number of particles
1000 atoms

# Specify the mass of each particle
mass 1 1.0

# Specify the initial positions and velocities of each particle
x y z vx vy vz
...

=========================================================================================================================

How to Simulate a Semi-Flexible Polymer in LAMMPS

In this blog post, I will walk you through the steps necessary to simulate a semi-flexible polymer in LAMMPS. I will start by explaining the basic concepts of polymer simulation, and then I will provide a step-by-step guide to setting up and running a simulation. Finally, I will discuss some of the results that you can expect to see from your simulation.

Q. What is a semi-flexible polymer?

A semi-flexible polymer is a long chain of molecules that is able to bend and twist, but which also has some resistance to these deformations. This resistance is due to the fact that the molecules in the polymer are connected to each other by chemical bonds. These bonds act like springs, and they resist being stretched or bent.

Q. Why simulate semi-flexible polymers?

Semi-flexible polymers are important in a wide variety of applications, including the design of new materials, the study of biological systems, and the development of new drugs. For example, semi-flexible polymers are used in the manufacture of everything from plastics to cosmetics. They are also found in many biological systems, such as single and double stranded DNA and proteins.

Q. How to simulate semi-flexible polymers in LAMMPS?

To simulate a semi-flexible polymer in LAMMPS,
1. You will need to create a system of particles that represent the individual molecules in the polymer.
2. You will also need to specify the force field that will be used to describe the interactions between the particles.
3. The force field will need to include terms that account for the bonded interactions between the molecules, as well as the non-bonded interactions between the molecules.

Q. How to create a system of particles that represent the individual molecules in the polymer

To create a system of particles that represent the individual molecules in the polymer, you will need to specify the following information:
1. The number of particles
2. The mass of each particle
3. The initial positions and velocities of each particle.
I use a simple python script to generate the data file. The script is available on my github (link).

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Define the variable list

The first step is to define the variables that will be used in the simulation. The variables are:
1. ParticleN: The number of particles in the system.
2. temp: The temperature of the system.
3. kappa: The spring constant of the bonds between the particles.
4. gamma: The friction coefficient of the system.
5. damp: The damping coefficient of the system.

# Variable List
variable ParticleN equal 64
variable temp equal 1.0
variable kappa equal 8.0
variable gamma equal 0.7 
variable damp equal 1.0/${gamma}
dimension 3
units lj
atom_style	molecular
#boundary p p p # periodic boundary
boundary s s s # fixed boundary
#boundary f f f # shrink-wrap boundary
comm_modify cutoff 4

#neighbor settings
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes

Read the particles' data

Next, read the data file (named as poly1_ParticleN) that contains the specifications of the particles.
After reading the data file, you can specify the bounding box for the simulation.

#=================================Load Initial Configuration=================================#
read_data poly1_${ParticleN}.input
#============================================================================================#

#Box Boundary
region box block -5000 5000 -5000 5000 -5000 5000

#Particle Properties 
mass * 1 # assign the mass of the monomer (for a homopolymer mass of the monomer is 1)) 
group polymer3d type 1 # assign the group polymer3d to all the particles

Assigning the initial velocity of the particles

The initial velocity of the particles is assigned using the command velocity.
The command velocity takes the following arguments:
1. The group of particles to which the velocity will be assigned.
2. The temperature of the system.
3. Random number seed.
4. The type of velocity distribution.

#========================================Initial Velocity====================================#
velocity all create ${temp} 1298371 dist gaussian
#============================================================================================#

Defining the interaction potential among the monomers

We use a short range repulsive Lennard-Jones (LJ) potential, \[U_{LJ}(r) = 4\epsilon \left[{\left(\frac{\sigma}{r}\right)}^{12}-{\left(\frac{\sigma} {r}\right)}^6\right]+\epsilon, \quad \textrm{for} \quad r\le 2^{1/6}\sigma\] \[= 0, \quad \textrm{otherwise}\] to model the excluded volume interaction between two beads separated at a distance r, where ϵ is the strength of the LJ potential.

The connectivity between two neighboring monomers is constructed using the finitely extensible nonlinear elastic (FENE) spring \[ U_{\mathrm{FENE}}(r_{ij})=-\frac{1}{2}k_FR_0^2\ln\left(1-r_{ij}^2/R_0^2\right). \] Here, \(r_{ij}=\left | \vec{r}_i - \vec{r}_j \right| \) is the distance between two consecutive monomer beads i and \(j=i\pm1\) at \(\vec{r}_i\) and \(\vec{r}_j\), \(k_F\) is the spring constant and \(R_0\) is the maximum allowed separation between two connected monomers.

An angle dependent three body interaction term is introduced between successive bonds which accounts for the chain stiffness \(\kappa\) \[U_{bend}(\theta_i) = \kappa\left(1-\cos \theta_i\right) \] and \(\theta_i\) is the angle between the bond vectors \(\vec{b}_{i-1} = \vec{r}_{i}-\vec{r}_{i-1}\) and \(\vec{b}_{i} = \vec{r}_{i+1}-\vec{r}_{i}\), respectively.

#===================================Potential Definition=====================================#
#Weeks-Chandler-Anderson
pair_style         lj/cut 1.12246204830937
pair_coeff         * * 1.0 1.0 1.12246204830937
pair_modify        shift yes

#FENE type bond
bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0
special_bonds fene

#Cosine type angle
#cosine shift angle: E = (U_min/2)*(1+cos(theta-theta_0))
variable K equal 2.0*${kappa}
angle_style cosine/shift
angle_coeff * ${K} 180.0
#=============================================================================================#

Remove the overlap of the beads (if any)

#remove overlap between polymer beads
minimize 1e-7 1e-7 1000 1000
reset_timestep 0

Langevin dynamics simulation

Once you have created your system and specified the force field, you can start the equilibration of the polymer. The simulation will proceed by solving the equations of motion for the particles. The equations of motion will take into account the forces acting on the particles, as well as the random forces that are due to the thermal motion of the molecules. The equilibration steps will be discarded and often called as pre-production run.

#####################################################
# Equilibration (Langevin dynamics at T)
#=======================================Set Integrator========================================#
fix integrator polymer3d nve gjf = yes
fix dynamics   polymer3d langevin ${temp} ${temp} ${damp} 252352

# specify timestep
timestep	0.01

# Thermodynamic Quantities Print 
thermo_style custom step temp  ke pe  
thermo_modify   norm yes flush yes
thermo    1000000

run 10000000 # Equilibriation Steps (will be discarded)
#=============================================================================================#

Post-equilibration simulation

# Reset the Integrator
reset_timestep 0
unfix integrator
unfix dynamics

#####################################################
# Post Equilibration (Langevin dynamics at T)
#=======================================Set Integrator========================================#
fix integrator polymer3d nve gjf = yes
fix dynamics   polymer3d langevin ${temp} ${temp} ${damp} 252352

# specify timestep
timestep	0.01

The post-equilibration simulation will continue for a specified number of steps. At each step, the positions and velocities of the particles will be updated. After the simulation has completed, you can analyze the results to obtain information about the properties of the polymer, such as its radius of gyration, its end-to-end distance, and its conformation.

################################################################################################
#                          On the Fly Statistical Quantities Calculations                      #
################################################################################################
# Radius of Gyration Output
compute 1 all gyration
variable Rgsq equal c_1*c_1
fix Rgave all ave/time 100 5 1000 c_1 v_Rgsq file Rg.dat 

# Potential Energy Output
compute 2 all pe
variable pesq equal c_2*c_2
fix Peave all ave/time 100 5 1000 c_2 v_pesq file PEavg.dat

# End to End distance
variable Rendsq equal (x[${ParticleN}]-x[1])^2+(y[${ParticleN}]-y[1])^2+(z[${ParticleN}]-z[1])^2
fix Rendave all ave/time 100 5 1000 v_Rendsq file Rend.dat

# Bond Length Calculation
compute 3 all bond/local dist 
compute Ave_BL all reduce ave c_3
variable BL equal c_Ave_BL

Output the thermodynamic quantities

# Thermodynamic Quantities Print 
thermo_style custom step temp v_BL 
thermo_modify   norm yes flush yes
thermo    10000

Output the trajectory

# Write Coordinates 
#dump 1 all atom 100 Running_Config.lammpstrj # (Movie Making)
dump 2 all dcd 5000 Running_Config.dcd 


run 10000000 # Post-Equilibriation Steps 

In this blog post, I have provided a brief overview of how to simulate a semi-flexible polymer in LAMMPS. Hope that this blog post has been helpful, and I encourage you look at the full code following the github link.

References

1. Thompson, A. P., Aktulga, H. M., Berger, R., Bolintineanu, D. S., Brown, W. M., Crozier, P. S., in ’t Veld, P. J., Kohlmeyer, A., Moore, S. G., Nguyen, T. D., Shan, R., Stevens, M. J., Tranchida, J., Trott, C., & Plimpton, S. J. (2022). LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Computer Physics Communications, 271, 108171.
2. Seth, S., & Bhattacharya, A. (2022). How capture affects polymer translocation in a solitary nanopore. The Journal of Chemical Physics, 156(24), 244902.