This is my documentation about LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator). Here I'm putting some very especific info for a very especific application. The official documentation cab be found in LAMMPS-OFF-DOC.
Input Script
- Initialization:
- Atom definition:
- Settings:
- Force field coefficients: pair_coeff, bond_coeff, angle_coeff, dihedral_coeff, improper_coeff
- Simulation parameters: neighbor, neigh_modify, group, timestep, reset_timestep, run_style, min_style, min_modify.
- Fixes to impose some boundary conditions, time integration, or diagnostic options.
- compute, compute_modify, and variable.
- Output options thermo, dump, and restart.
- Run simulation:
Commands: units, dimension, processors, boundary, atom_style, atom_modify
Create atoms on a lattice using lattice, region, create_box, create_atoms. Atoms cab be defined also froam a data or restart file.
Once atoms are defined , we need to set:
Run simulation using run. Energy minimization with minimize command.
Simple Script
# 2d Born potential gas
# U = A exp(\frac{\sigma-r}{\rho} - \frac{C}{r^6} + \frac{D}{r^8}
units lj
dimension 2
boundary p p p
atom_style atomic
# create geometry
lattice sq 0.8
region box block 0 35 0 35 -0.25 0.25
create_box 1 box
create_atoms 1 box
mass 1 1.0
# velocity all create temperature random_seed
velocity all create 1.0 200
velocity all zero linear
# Potential
pair_style born 15.0
pair_coeff * * 1.0 1.0 0.0 0.0 0.0
# fix
fix 1 all nve
# Run
timestep 0.001
thermo 1000
thermo_style custom step etotal ke pe temp press
dump 1 all custom 1000 dump.vel id vx vy vz
compute myRDF all rdf 50
# fix ID group-ID ave/time Nevery Nrepeat Nfreq
fix 2 all ave/time 1 100 1000 c_myRDF[1] c_myRDF[2] file tmp.rdf mode vector
run 10000
- Initialization:
- Atom definition:
- Settings:
- Run:
Units are Lennard-Jones (defalult value): all quantities are unitless. Mass, sigma, epsilon, and Boltzmann contant are set to 1. Other options are real, metal, si, cgs, electron, micro, and nano.
We are defined a two dimensional simulation with the command: dimension 2.
For our simulation we are using periodic boundary condition on all coordinates (defalult values). Other options are f: non-periodic and fixed, s: non-periodic and shrink-wrapped, and m: non-periodic and shrink-wrapped witn minimum value.
Atom style is atomic (defalult value), which is used coarse-grain liquids, solids and metals.
lattice sq 0.8: This is a 2D lattice with reduced density equals to 0.8 (arbitarily choosen), where $$ \rho^{*}=\frac{N}{V^{*}}=\frac{N}{{r^{*}}^d} $$ with \(N\) the number of basis molucules, \(V^{*}\) is the reduce volume, and \(d\) is the dimension.
region box block 0 35 0 35 -0.25 0.25: this simulation will create 1225 atoms, separated by 1 lj-unit, so in one of the directions we will have \(\sqrt{1225}=35\) atoms separated one cell unit. Then this command will create a box with coordinates in x from 0 to 35, in y from 0 to 35 and from -0.25 to 0.25 in z. The input box is just an ID asssigned by the user, and block is a style of region (others options are cone, cylinder, etc).
create_box 1 box creates simulation box for N types of atoms (1 in this case) for the region with ID called "box".
create_atoms 1 box creates N types of atoms (1 in this case) in the region with ID called "box". The simulation box must already exists.
mass 1 1.0 Set the mass for the atom type N (in this case mass=1.0 for atom type 1).
velocity all create 1.0 200 creates ensamble of velocities for 'all' atoms using random number generator (with seed=200) at a specific temperature (T*=1.0).
velocity all zero linear set the total linear momentum of the atomic system to zero.
pair_style born 15.0
pair_coeff * * 1.0 1.0 0.0 0.0 0.0
Here we will use the Born-Mayer-Huggins potential, with cutoff radius of \(r_c=15.0\) lj-units. The potential is given by $$ E = A \exp\left( \frac{\sigma - r}{\rho} \right) - \frac{C}{r^6} + \frac{D}{r^8}, ~~r < r_c $$ For now, we just need the exponential part of the potential so we set its parameters with pair_coeff: first two entries specified the atoms with this interaction (in this case the interaction is for all atoms, so we used * *), then we specified the numerical values of the constants A=1.0, \(\rho=1.0\), \( \sigma=0.0 \), C=0.0, and D=0.0
fix 1 all nve The fix 'nve' is applied to a group of atoms ('all' in this case). nve performs a constant NVE integration, this creates a system trajectory consistent with the microcanonical ensamble.
timestep 0.001 Set timestep
thermo 1000
thermo_style custom step etotal ke pe temp press Output themodynamics info every N steps (1000 in this case-- If 0 just gives info in the initial and final timestep), In thermo_style we specified the info we want to output (custom). In this case we are getting the number of timestep, total energy, kinetic and potential enegies, temperature and pressure.
dump 1 all custom 1000 dump.vel id vx vy vz Dump a snapshot of atoms N (1 in this example) quantities (velocities components in this case). The snapshots are stored in the file 'dump.vel'
compute myRDF all rdf 50
Compute the Radial Distribution Function for all particles with N=50 bins
fix 2 all ave/time 10 100 1000 c_myRDF[1] c_myRDF[2] file tmp.rdf mode vector
The averaged RDF is generated on timesteps multiple of Nfreq. The average is calculated over Nrepeat quantities, computed using samples taken every Nevery.
Nfreq must be a multiple of Nevery. Nrepeat*Nevery can not exceed Nfreq. If Nrepeat=1 and Nfreq=100, there are not time averaging, values are generated on timesteps 100, 200, etc
run 10000 Run the simulation for 10,000 timesteps.