This page provides general information about the functionnement of the open source software LAMMPS, which according to its own website is a classical molecular dynamics (MD) code that models ensembles of particles in a liquid, solid, or gaseous state.
Warning : this is not an exaustive list of the possibilities offered by LAMMPS, for that see the LAMMPS documentation.
Credit : This page has been inspired by this github post by Ryan Kingsbury.
Input files
Any lammps simulation requires at least one input script containing commands. In addition, data files containing information about the system topology can be supplied, as well as parameters file containing information about the force field, such as masses and pair coefficients.
The input script
Input scripts are text files that are read by LAMMPS one line at a time. Input scripts contain commands (see below), and each command induces LAMMPS to take action. Line starting with # are ignored by LAMMPS. This is an example of a minimal script extracted (with some simplifications) from this tutorial:
# Minimal LAMMPS input script with comments (this line is ignored)
pair_style lj/cut 2.5 # use Lennards Jones potential with a 2.5 (unitless) cutoff
region myreg block -10 10 -10 10 -10 10 # define a cubic region
create_box 1 myreg # set the simulation box according to the region "myreg"
create_atoms 1 random 100 54142 myreg # create 100 atoms of type 1 in the region "myreg"
mass 1 2 # attribute a mass of 2 (unitless) to atoms of type 1
pair_coeff 1 1 1.5 0.7 # attribute Lennards Jones coefficients (epsilon = 1.5 and sigma = 0.7) to atoms of type 1
dump mydmp all atom 10 dump.lammpstrj # dump the atoms position in a file every 10 timestep
minimize 1.0e-4 1.0e-6 1000 10000 # perform an energy minimisationWhat is a fix ? A fix is a command used to apply an operation on a group of atoms. For instance, thermostating can be applied using fix Berendsen/temp or fix temp/rescale, integration of atom positions can be done using fix nve, or an additional force can be applied to the atoms using fix gravity.
Can I create a loop ? Yes, ‘for loop’ can be created usign the jump command, and ‘while loop’ can be created combining the jump command with the if statement, see this post.
System creation
Typically, a simulation requires the definition of a computational box (orthogonal or tilded), as well as a list of atoms and their initial positions. These information can either be written in a data file before the simulation is started, or be specified within the input script.
Option 1: from a data file
LAMMPS can read information from a data file using the read_data command. The data file usually contains information such as the number of atoms and their initial positions, the box dimensions, etc. Such data file must follow a specific formating, see here. Here is an example of a data file for a simple system with 1 water molecule:
15 atoms
2 atom types
2 bonds
1 bond types
1 angles
1 angle types
-5.00 5.00 xlo xhi
-5.00 5.00 ylo yhi
-5.00 5.00 zlo zhi
Masses
1 15.9994
2 1.008
Atoms
3 1 2 0.527 -0.43213494211794057 4.490986821056591 -2.00832253672005
1 1 1 -1.054 -1.0103199635859141 3.9862944633721575 -1.436277429248943
2 1 2 0.527 -1.819251731628718 3.888533765608223 -1.938580267494636
Bonds
1 1 1 2
2 1 1 3
Angles
1 1 2 1 3Some external softwares allow you to create such data file, for example VMD + topotool. Data files can also be created with a simple script, find some Python and Octave/Matlab examples in my Github repository.
Option 2: LAMMPS build-in commands
Another way to create a simulation box and place atoms is to use the create_box and create_atoms commands of LAMMPS. For instance, to generate a box of 10x10x10, type in the input script:
region box block -5 5 -5 5 -5 5To create a simulation box with the space for 2 atom type, 1 bond type, and 1 angle type, write in the input script:
create_box 2 box &
bond/types 1 &
angle/types 1 &
extra/bond/per/atom 2 &
extra/angle/per/atom 1 &
extra/special/per/atom 2And to place a single water molecule, type:
molecule h2omol H2OTip4p.txt
create_atoms 0 random 1 456415 NULL mol h2omol 454756Here the molecule template H2OTip4p.txt contains information about the molecule (i.e. which atoms are bonded toghether, etc.).
Outputs
There are several ways to extract information from a simulation:
Thermodynamic info
Thermodynamic info (such as temperature, energy, pressure) can be printed in the terminal and in the log file by using the thermo command. To print thermodynamic information every 100 timestep, type:
thermo 100The thermo_modify and thermo_custom commands can be used to change the type of info printed by LAMMPS, for instance, to print the temperature of the group groupH2O only, as well as the overall kinetic energy and volume, use:
group groupH2O type 1 2
compute tempH2O groupH2O temp
thermo_modity temp tempH2O
thermo_style custom step temp ke volTrajectory
Trajectory information, i.e. the position of atoms over time, can be printed in a file using the dump command:
dump myDump all atom 100 dump.lammpstrjHere the positions of all the atoms are printed every 100 timesteps in the dump.lammpstrj file, which you can then open using VMD or Ovito. ‘lammpstrj’ is the default ascii LAMMPS dump format. Compressed trajectory formats such as xtc can also be used.
There are plenty of measurements that can be performed from the trajecory files after your simulation is over. These measurements can be done with external tools such as MDAnalysis or MAICoS (for confined geometries).
Data
Data can be printed in a text file as a function of the simulation timestep using fix ave/time. The printed information can for instance be the result of a variable, a compute, or even a fix. For instance this series of commands extracts the the box size along x (lx), convert it in meter, and print it in a text file every 1000 timesteps:
variable AngstromtoMeter equal 1e-10
variable mylx equal lx*${AngstromtoMeter}
fix myat1 all ave/time 1000 1 1000 v_mylx file myly_in_meter.datData can also be printed as a function of space, for instance to extract the 1D density profile along z, type:
compute cc1 water chunk/atom bin/1d z 0.0 1.0
fix myac1 water ave/chunk 100 1000 100000 cc1 density/number norm sample file density_water.dat
