LAMMPS: An interatomic potential energy code

We will be using the LAMMPS simulation software. LAMMPS allows performing a variety of simulation types on periodic and isolated structures. Further information (including an extensive online manual) can be found on the LAMMPS website:

http://lammps.sandia.gov

Manual: http://lammps.sandia.gov/doc/Manual.html

LAMMPS is free for academics and provides several interatomic potential models for a broad range of systems.

Ask When You get Stuck

If you run into questions, email Alex Urban (au2229@columbia.edu) or TA Xinhao Li (xl2778@columbia.edu), post on Slack, or stop by during office hours.

Note: all questions should be submitted at least 24 hours before the project is due in order to guarantee a response.

This Project uses nanoHUB.org

To run the simulations for project 1, we will make use of the open nanoHUB infrastructure. Please register (for free) on https://nanohub.org to be able to access the simulation tools. We will use the LAMMPS tool to carry out all simulations for project 1.

Problem 1

Read over problem 1 in the assignment.

For problem 1, we provide LAMMPS input files to get you started. You will have to adjust these files yourself for the subsequent problems. Generally, LAMMPS requires at least one  input  file  defining  the  structure  (atomic  positions,  lattice),  energy  model,  and calculation details.

LAMMPS input files

Let us take a look at the LAMMPS input files for problem 1a. The file ‘lmp-in.1a- single’ specifies the parameter for a single point calculation (energy evaluation) of Pt in the FCC structure.

The first lines in the input file take care of the general set-up:

The next block in the input file creates a face-centered cubic (FCC) Pt structure:

0 orient y 0 1 0 orient z 0 0 1 0 1  0 1  0 1 units lattice

create_box

create_atoms

replicate

create an FCC lattice unit cell with a lattice parameter of 4.0 Ǻ and lattice vectors that are oriented parallel to the Cartesian directions.

create a region named ‘simulation_box’ that covers the entire unit cell. set the region ‘simulation_box’ to be used as simulation box.

decorate the lattice sites with atoms in the entire simulation box.

construct a 1×1×1 supercell. The size of the supercell will have to be adjusted

for the latter problems.

The next line simply defines the mass of an Pt atom in atomic mass units:

11.   mass 1 195.084

The following  block  contains  the  definition  of the  energy  model,  i.e.,  the  type  and parameters of the interaction potential.

12.   pair_style lj/cut 10.0

13.   pair_coeff 1 1 0.200 2.540

Line 12:  type of the interaction potential. For problem 1a this is a Lennard-Jones potential with a cutoff distance of 10 Å

Line 13:  parameters  of  the  Lennard-Jones  potential  for  the  interaction  of  atoms  of

species 1 (Pt) with each other: ε = 0.200 eV and σ = 2.540 Å        Note that LAMMPS uses the following form of the Lennard-Jones potential:

VLJ (r) = 4 e  )*  , !2  − *  , 6 .

The next lines in the input file specify details of the simulation and actually carry it out. In ‘lmp-in.1a-single’ this is simply an energy evaluation without any modification of the structure:

14.   fix 1 all nve

15.   thermo_style custom step pe lx ly lz

16.   run 0

Line 14:  in LAMMPS terminology, a ‘fix’ is an operation that is done at every step of the simulation. The ‘nve’ fix takes care of evaluating the potential energy.

Line 15:  the ‘thermo_style’ defines the quantities that should be reported in the output file. Here, we request the potential energy (pe) and the length of the simulation cell in the three dimensions (lx, ly, and lz)

Line 16: run zero dynamics steps. This simply evaluates the energy of the initial structure. In later projects, you will learn how to run molecular dynamics simulations using

LAMMPS.

The final two lines of the input file are only for your convenience. They just print out the final potential energy on a separate line in the output file.

Relaxations

The second input file for problem 1a, ‘lmp-in.1a-relax’, uses the same structure and energy model definitions, but specifies a different kind of simulation. Here, LAMMPS is used to perform a structure relaxation, i.e., to optimize the atomic positions and the cell parameters.

The lines with the simulation parameters now read:

17.   fix 1 all box/relax iso 0.0 vmax 0.001

18.   thermo 1

19.   thermo_style custom step pe lx ly lz

20.   minimize 1e-15 1e-15 5000 10000

Line 17:  define a ‘box/relax’ fix that takes care of the relaxation.

Line 18:  this  line  only  affects  the  output  file.  ‘thermo  1’  makes  LAMMPS  write  out thermodynamic quantities after every 1 simulation steps (i.e., at each step).

Line 19:  is identical to line 15 above.

Line 20:  carry out an energy and force minimization. The convergence criteria for both energy and forces are set to 10- 15 (relative) and 10- 15 eV/Ǻ, respectively. We limit the total number of optimization steps to 5,000 and the total number of force evaluations to 10,000.

EAM potential

The two input files discussed so far both define a Lennard-Jones potential as energy model. For problem 1b, you will need equivalent input files for an embedded atom model (EAM) potential.

As you can see in ‘lmp-in.1b-single’ and ‘lmp-in.1b-relax’, only the lines related to the potential definition have to be updated and now read:

21.   pair_style eam

22.   pair_coeff * * Pt-Adams1989.eam Pt

Line 21:  use an EAM alloy-type potential.

Line 22:  load EAM parameters for Pt from file ‘Pt-Adams1989 .eam’ .

Of course, the EAM potential file ‘Pt-Adams1989.eam’ has to be available for LAMMPS to use it. You will have to upload it to nanoHUB along with the general LAMMPS input file.

Running LAMMPS

To run LAMMPS using any of the provided input files, first open the LAMMPS tool on the nanoHUB website. Next, upload the input file via right click on the input file field. For the EAM potential simulations, also upload the EAM potential file on the ‘Simulation’ tab (activate ‘Upload additional files’) and make sure that the file name is correct. Finally, click on the ‘Simulate’ button to start the simulation. All simulations for project 1 should only take a few seconds, so if the simulation takes more than a minute ‘Abort’ it and check your input file.

LAMMPS output files

Once  the  simulation  is  done,  the  nanoHUB  LAMMPS  tool  will  let  you  inspect  the generated output files. The most important output (for project 1) is the stdout file which contains the output that would usually be printed to the screen during the simulation.

Select this file in the ‘Result’ menu.

Skip over the stdout file until the line:

Step PotEng Lx Ly Lz

The subsequent lines contain the ‘thermo’ output that we requested in the input file, i.e., the potential energy (PotEng) and the simulation cell dimensions (Lx Ly Lz).

At the very end of the stdout file you will find the final total (potential) energy.

Problem 2

Read over problem 2 in the assignment.

For this problem, you will have to modify the input files provided for problem 1.

To determine the converged vacancy formation energy you will have to create supercells of the FCC unit cell. This is simply done by changing the line

replicate 1 1 1

in the input file.  For example, to create a 3×3×3 supercell, the command would be ‘replicate 3 3 3’ .

One way to generate a vacancy is by removing the first atom from the structure. This can be accomplished in LAMMPS with the following additional lines in the input file (to be inserted after the supercell creation):

23.   group vacancy id <= 1

24.   delete_atoms group vacancy

Line 23:  create an atom group named ‘vacancy’ containing all atoms until index 1. This means only a single atom will be assigned to this group.

Line 24:  delete all atoms within the group ‘vacancy’ .

As for problem 1, don’t forget to upload the EAM potential file for those simulations that use it.

Problem 3

Read over problem 3 in the assignment.

Assume a simulation box full of atoms has been created in the following way:

# Create atoms

To build a slab model to calculate surface energies, add one line in the LAMMPS input file to create a vacuum region:

change_box all z delta 0.0 0.5 units lattice

The above line shifts the simulation box along the z direction by 0.0 lattices for the lower bound and by 0.5 lattices for the upper bound. This shift creates a vacuum region of 0.5 lattices, i.e., 2 Ǻ in the example.

Calculate the energy with and without the vacuum. Use the properly normalized energies to calculate the surface energy.

The structure relaxation can be omitted.

FAQ for project 1

What dimensions of supercell should I use?

Remember periodic boundary conditions. This is the most important.

For the surface energy calculations, the cell will be different in the z direction. Remember though, because of periodic boundary conditions, there is no need to enlarge the cell in x or y direction. So 1×1×5 or 1×1×10 are suitable supercells. You could also build a 2×2×10 cell, but this would just increase the compute time and would not improve the result over the 1×1×10 cell. There are two convergence issues to think about, the slab (the part with the atoms) thickness, and the vacuum thickness.

Note, this cell dimension choice does not apply to all problems! Before you choose your supercell, you have to think about the problem you are working on. For instance, it the situation is different for the vacancy problem! Think about why this is.

Lastly, you cannot know for sure beforehand what an appropriate supercell size is. You must always test for convergence. For example, if you look for the vacancy formation energy in a 2×2×2 supercell, you must test a larger supercell to make sure that your supercell is big enough. And the required supercell size can differ from potential to potential.

In   the   vacancy   problem,   is   the   answer   I   get   only   applicable   to   the   vacancy concentration  that  I  have  created?  I  do  not  understand  the  supercell  convergence issues for the vacancy problem .

Remember the objective is not to simulate the real material, but to compute the properties  that are relevant for the real material. For example, say one has a metal with vacancy  concentration 10-6 , which is quite typical for real metals. Does this mean that one should  study the vacancies in such a metal with a supercell with 1 million atoms and one vacancy? Of course not! At such low concentration the vacancies are isolated. (That means that  they do not interact. They do not "see" each other). So any vacancy that behaves isolated  is the same, whether it is at concentration 10-8 or 10-5 . What we need to do is find the  smallest supercell where it is still a reasonable approximation that the vacancies in the  different images of the cell do not interact. One way is to make the cell systematically  bigger and see if the energy converges. When the vacancy formation energy does not  change much anymore as the cell gets bigger, one has found the formation energy for an  isolated vacancy. In practice this occurs at much higher concentrations than 10-6 (and  hence for smaller supercells).