Write Modular code in LAMMPS.

4 min readFeb 6, 2023

Lifeless Introduction

LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) is an open-source software for molecular dynamics(MD) simulations. It provides a platform for simulating the physical movements and interactions of atoms and molecules in materials, enabling the prediction of their properties and behavior. LAMMPS is widely used for simulating materials science, condensed matter physics, biophysics, and many other areas in which molecular-level understanding is crucial. The software is designed to be highly scalable and capable of running on large parallel computing platforms, making it well-suited for simulating complex, large-scale systems.

Ok, Modular code, but why?

Writing Modular code can help to rerun a specific part of the simulation or easily restart a failed simulation without starting it from the beginning

I think the reason I saw most people don’t write modular code in LAMMPS is that LAMMPS documentation provides examples of codes which is not modular. But we can follow the same approach we use in GROMACS (For those who are not familiar with GROMACS, it is another MD simulation package that is mostly used for biomolecular simulations) to write LAMPPS code which will help us to easily restart or rerun a specific part of our entire simulation.

Benefits of Writing Modular Code in LAMMPS

Reusability: Modular code is easier to reuse because it is organized into separate units that can be combined to create a complete simulation. This makes it easier to reuse code for multiple simulations and reduces the time spent writing new code.
Maintainability: Modular code is easier to maintain because changes to one unit will not affect the rest of the code. This makes it easier to update the code as needed and reduces the risk of introducing bugs.
Readability: Modular code is easier to read because it is organized into separate units, each with a well-defined purpose. This makes it easier for others to understand the code and collaborate on its development.

I am convinced. How to do it?

You are here, so I assume you already know the basic structure of the LAMMPS code. If you are not, please refer to the LAMMPS beginner guide to get started.

LAMMPS BEGINER HELP

So our code will contain 4 files. The input script. The initial file, the settings file, and the data file.

The Input Script

This is the main file you will pass to LAMMPS, and for each step, i.e., energy minimization, energy equilibration, and deformation, you will have different input scripts. One sample input script for Energy minimization can be like this:

#### Minimization

################################ SYSTEM INIT #################################
include initial.in
read_data initial.data
include settings.in

################################## OUTPUT ########################################## 
thermo_style    custom step pe etotal vol epair ebond eangle
thermo          20  

################################# Minimization ################################# 
minimize 1.0e-8 1.0e-10 10000 300000

########################################################################################
write_data  system_after_min.data
########################################################################################

see how we are loading the initial.in and settings.in file in our main initial file. We are describing the atom structure in initial.data file. The structure of the initial.data file can be different based on the force field you used.

The most important part is to see how we save the data file after minimizing the structure. This file can be used in the next step of energy equilibration by replacing the input file while keeping the initial and settings file unchanged.

Here is an example of how we are starting the simulation from the output of the previous step by reading the structure using read_data system_after_min.data command.

####RUN INPUT FILE
################################ SYSTEM INIT #################################
boundary p p p
include initial.in
read_data system_after_min.data
include settings.in
timestep      0.01
velocity all create 1.0 12345

########################################################################################
group lipid type 1 2
group bolalipid type 3 4 5
group network type 6 7 8
group spectrin type  7 8

################################## OUTPUT ##########################################
dump 1 all custom 2000 traj_npt.lammpstrj id mol type x y z ix iy iz
thermo_style  custom step pe etotal vol epair ebond eangle
thermo        200  # time interval for printing out "thermo" data
############################### Equilibration ################################# 
fix fxlan all langevin  1.1  1.1  1.0   48279
fix fxnve all nph iso 0.0 0.0 10
########################################################################################
run 100000
write_data system_after_equil.data

The Initial File

This file will contain the initialization of the MD simulation required for running any MD simulation. See an example:

# Initializing the system

units           lj
atom_style      full
bond_style      hybrid fene 
angle_style     hybrid harmonic  
pair_style      hybrid lj/cos2 1.12246204831 lj/wca 1.12246204831
special_bonds   lj 0.0 1.0 1.0
neigh_modify    every 1 delay 1
neighbor        1.6 bin

The Settings File

This file contains the description of the system and force field coefficient. Here is how I do it:

#PAIR COEFFICIENTS
pair_coeff  1 1 lj/wca 1 0.95
pair_coeff  1 2 lj/wca 1 0.95
pair_coeff  2 2 lj/cos2 1 1 1.6    

#BOND COEFFICIENTS
bond_coeff  1    fene      30.0 1.5 1.0 1.0       #Lipid

#ANGLE COEFFICIENTS
angle_coeff 1  harmonic  10.0 180          #Lipid

I am not describing the data file here as it can be very specific to the forcefield.

Disclaimer:

This input script is solely for demonstrating how LAMMPS code can be written in a modular format. The parameters used here are not the parameters I am suggesting you should use. Please follow the literature and your judgment to select parameters for your specific application.

Conclusion

Thank you for taking the time to reach out all the way here. If you want to thank me for my help, you can do so by hitting the follow button or citing my paper.

Khayrul Islam, Sourav Saha & A. K. M. Masud (2020) Molecular dynamics simulation of the mechanical properties of CNT-polyoxymethylene composite with a reactive forcefield, Molecular Simulation, 46:5, 380–387, DOI: 10.1080/08927022.2020.1711904