LAMMPS Nanoindentation Tutorial

Directory


Disclaimer

Despite the fact that this tutorial may make me sound like an authority (although I try to remember to put comments in places where I am unsure), I am far from an authority on LAMMPS. I didn't write it, and I am not affiliated with its authors or Sandia National Lab in any way. I am only a user. I think it works very well, and I like using it, but there may be other codes that are better for your purposes out there.

Generating LAMMPS Input

LAMMPS has a very well written manual, in general, but I think there are some issues that need to be clarified espically with regards to performing Molecular Mechanics calculations. In this document, I use a monospaced typewriter typeface to explicitely reference commands, folders, and files.

Preliminaries

Units In LAMMPS, there are two meanings of units: the unit system used when calculating physical quantities, and the units in which to specify distances (either in fractions of a unit cell, or in absolute units which are always in Å). LAMMPS is inconsistent with the default units for the specification of distances, espically with the construction and fix commands. Regardless of what the manual says, always specify the units on the commands that take units, even if the command has default units.

A critical feature that will be put to extensive use in our model is the fix function. Fixes are used to apply forces, velocities, or other constraints on a specific group of atoms. The fix command, however, only acts on groups of atoms that have been defined by the group command. See the LAMMPS manual to learn about defining groups as I use them in the nanoindentation input files.

Writing Nanoindentation Input Files

This section describes how to write LAMMPS input files for use with the procedures and scripts described in the discussion of visualization using VMD below. My procedure for modeling nanoindentation consists of two steps:
  1. Create and relax initial configuration
  2. Perform indentation beginning from relaxed configuration
In order to save time, we can use LAMMPS to create a restart file for the relaxed initial configuration that can be used repeatedly as a starting point for simulations on structures of the same geometry. The example which will be continually discussed on this website will be a hemispherical Nickel nanodot with a diameter of 100Å, fixed to a ridgid substrate.

Establish the Initial Configuration

The first step to model the indentation of a nanodot is to create the initial geometry. In the restart_100A.lmp input file, the geometry of a hemispherical nanodot with a diameter of 100Å. The next step is to find the initial relaxed configuration of the nanodot.

The input file for the actual simulation is actually much simplier to write. In essence, you need to first load the relaxed nanodot configuration which will be found in a file (if using my input file) called: restart_100A.###, where the pound signs indicate a number that corresponds to the number of steps that it takes to find the final configuration.

Perform the Indentation

The input file for the actual simulation is actually much simplier to write than the geometry file. In essence, you need to first load the relaxed nanodot configuration which will be found in a file (if using my restart_100A.lmp input file) called: restart_100A.###, where the pound signs indicate a number that corresponds to the number of steps that it takes to find the final configuration. I have provided the indent_100A.lmp input file as an example.

Each indentation step consists of three lines specifying the current distance of the indenter, a line telling LAMMPS to include the energy of interaction between the indenter and the structure (otherwise there would be no energy change in the system) and the minimizer command itself. Since this command must be repeated at each distance step as the indenter moves, I have devised the stepgen_inout.awk script to automatically generate the necessary indentation and withdraw steps which can be pasted into your LAMMPS input file. The use of the script is self explanatory.

Running Simulations

Molecular Mechanics

The default build of LAMMPS (as of 22April08) uses low tolerances for minimizing energy and force. The system energy often converges very quickly and is nearly always the stopping condition of the minimization step. The 22April08 version of LAMMPS allows for independent specification of both force and energy convergence criteria. LAMMPS will terminate the convergence if any of three conditions are met: energy convergence, force convergence, or a failure of the linesearch algorithm. It is important to understand the exact convergence criteria used:
Energy
The energy convergence criteria is based on having the average of the energy of the last two iterations times the specified convergence criteria being less than the change in energy between the last two iterations.
Force
The force is considered converged when the 2nd-norm of the force reaches the specified tolerance. Note that this quantity is dependent on the number of atoms in the system.

Best Practices

After a lot of experimentation with multimillion atom systems, you should set the energy convergence near machine zero, which is approximately 2.3 × 10-16 for Sapphire, since the energy criteria is much easier to meet than the force criteria. In addition, the energy criteria does not guarantee that even a local energy minima is achieved. For the force criteria, my experience indicates that the best practice is to actually modify the code to place the convergence criteria on the maximum component of the force, rather than the 2nd norm. With this condition, at least, we have a consistently size independent measure of convergence. For the simulations that I was conducting, I modified the Conjugate Gradient minimizer that was included with LAMMPS. My version of min_cg.cpp is here and is modified in lines 303-327.

Submitting Jobs

To avoid squabbling over individual nodes, we'll let the computer distribute the load for us. We are using ROCKS Linux cluster management software to accomplish this goal. I have developed a simple script where you can submit a lammps input script with my bsub-lammps-sapphire.sh script.
Note that in order to get LAMMPS to work, you need to have both the input file and a copy of your LAMMPS executible in the same directory. You must also execute your bsub job submission script in that same directory. If you modify the provided input script, LAMMPS (or bsub) may have trouble finding input unless you are carefull. Don't panic, first follow my troubleshooting steps below.

Troubleshooting:
  • Make sure that you have both your LAMMPS executible and input file in the same directory
  • Execute your bsub job submission script in the same directory in which your executible and input file reside
  • Specify the full path to your executible, e.g. $PWD/lmp_sapphire, but not to your input file (I don't know why)

Post-processing

Output Parsing

get_force.sh script to extract the position, and force on the indenter and the potential energy of the system into a single file.

Visualization

Visualizing the Centrosymmetry Parameter in VMD

Specifically, I am interested in visualizing the centrosymmetry parameter of Kletchner, et al. which is a useful way to visualize the dislocations in a crystal with a centrosymmetric point group.

In order to perform the visualization, I have quite a few scripts and settings the must be executed. I first include the command,

dump any_name all custom 20000 any_file.csp tag type x y z c_CSP
in my LAMMPS input file to create an output file that reports atom labels, coordinates, type, position and the centrosymmetry parameter of each atom.

Then, I employ the  script to convert the raw output of LAMMPS into the Protein Databank format (PDB) that is easily used with VMD.

A link off of the VMD website pointed me to a script that can color atoms based on a scalar parameter. The pdbofactor.tcl script loads a scalar present in the bfactor field of a PDB file into the user field of VMD. To load the parameter into VMD's user field, which is one of the ways to color atoms, you must first load the script into VMD before executing it by issuing the

source pdbofactor.tcl
command in the TCL console (found in the VMD 'Tools' menu). To avoid doing this each time you start VMD, you can add it to your .vmdrc file which is automatically run when you start VMD. This file analogous to your shell's profile file. The .vmdrc file should also be located in your home directory (if not, make one). All you need to do is add the 'source' command to the file.

Now, invoke the pdbofactor.tcl script from the TCL console with its argument being the PDB file you want to visualize. In the case of the present example:

pdbofactor indent_100A.pdb

Finally select user in the trajectory  submenu of the color atoms by ... drop-down box as shown below. The atoms in the structure should now be colored. To visualize the dislocations internally (using the indent_100A.pdb file), simply recreate the settings in the screenshots below.

You should get:

Remember, the settings shown below correspond the results of running the indent_100A.lmp which is made of Nickel. The actual value of the centrosymmetry parameter that corresponds to surfaces, partials, and stacking-faults is dependent on the lattice constant and structure. For the case of Nickel, it was determined that values ranging from 0.1 to 3 correspond to partial dislocations (dark red), and values between 3 and 7 correspond to stacking faults (light-red to nearly white). Atoms having CSP values greater than 10 are either at or near the surface of the nanodot. CSP values in the range of 10 to 13 (white) tend to occur when atoms at on a step on the surface of the nanodot, and those near 20 are surface atoms in a plane (dark blue). To shear off the surface atoms and set the color scheme to match the example given, use the Selected Atoms settings:

user > 0.1 and user < 20 and ((x-80)^2 + (y-80)^2 + (z-80)^2) < 47.5^2
For Best Results: The key to peering inside the structure is to eliminate the surface atoms. For Nickel, subtracting 2.5Å from the radius is necessary to completely eliminate the atoms on the surface (in the example above, the radius is 50Å). The center of the dot happens to be at x=y=z=80Å in this case.

Interpreting CSP Values and Setting Color Range

You may be wondering how I knew what CSP values to use as the range to display as the origonal paper by Kletchner, et al. only gives the range and type of dislocations for Gold. In theory, the proper CSP values could be calculated, but there is always going to be a distribution about the ideal values for stacking faults and partial dislocations.


 This data is for a Nickel dot with lattice constant of 2.54Å at 0K.

I have concluded that the best way to avoid this difficulty is to take an arbitrary frame from your PDB file of a simulation and create a histogram of the population of CSP values as shown below. Observe that there are four significant peaks from 0 to 20. The peak near zero has a very large number of atoms (number is too large to be shown on this plot) which corresponds to an insignificant degree of disorder from a perfect lattice in the interior of the dot. The next peak, will correspond to partial dislocations. The third peak from the left corresponds to atoms in volved in stacking faults. Atoms with CSPs lying in the the fourth peak or greater correspond to atoms which are at or near the surface of the dot.

You can verify my conclusions yourself by using VMD's query function to obtain an atom number and look up it's CSP value in the PDB file.

Movies

Error Messages
P4_error
This is a segmentation fault, or out-of-memory error. There are two sources of this error. It is most often caused by trying to write data outside the bounds allocated to an array. In LAMMPS, this is likely due to trying to use the create atoms command outside of the simulation domain box, as defined using: create_box region 1 domain. This error is also known to occur when a certain individual (who will remain unnamed) runs memory intensive simulations on the head-node.
127
Most often due to specifying the incorrect input filename in the bsub job submission script.
ċ
66stepcsp.txt
(139k)
Christopher O'Brien,
May 22, 2010, 2:14 PM
ċ
VMD_100A_script.tcl
(14k)
Christopher O'Brien,
May 22, 2010, 2:12 PM
ċ
bsub-lammps-sapphire.sh
(1k)
Christopher O'Brien,
May 22, 2010, 2:14 PM
ċ
csp2pdb.awk
(1k)
Christopher O'Brien,
May 22, 2010, 2:01 PM
ċ
data2xyz.awk
(1k)
Christopher O'Brien,
May 22, 2010, 2:01 PM
ċ
force_indent_100A.txt
(3k)
Christopher O'Brien,
May 22, 2010, 2:14 PM
ċ
getforce.sh
(1k)
Christopher O'Brien,
May 22, 2010, 2:01 PM
č
inden_100A.mpg
(1197k)
Christopher O'Brien,
May 22, 2010, 2:02 PM
ċ
indent_100A.lmp
(16k)
Christopher O'Brien,
May 22, 2010, 2:14 PM
ċ
indent_100A.log
(129k)
Christopher O'Brien,
May 22, 2010, 2:13 PM
ċ
restart_100A.lmp
(2k)
Christopher O'Brien,
May 22, 2010, 2:12 PM
ċ
stepgen-inout.awk
(1k)
Christopher O'Brien,
May 22, 2010, 2:12 PM
Comments