CAVS / Capabilities / Materials Modeling / ICME / Code Repository / LAMMPS Tutorial 2

Integrated Computational Materials Engineering (ICME)

LAMMPS Tutorial 2

Abstract

This is part two of a tutorial for running LAMMPS simulation on a Windows machine. For this example, the molecular dynamics simulation calculates the equilibrium lattice constant and the corresponding cohesive energy for aluminum. This tutorial will introduce the use of variables via command line and MATLAB for running LAMMPS.

Author(s): Mark A. Tschopp

Complete Tutorial One

If you have not done so already, complete the first tutorial available here.

Modify the input file

This input script was run using the Aug 2015 version of LAMMPS. Changes in some commands may require revision of the input script. In the below input script, the initial lattice constant (4) was replaced with the variable latconst (${latconst}). Either modify the initial script 'calc_fcc.in' or copy the text below and paste it into a new text file. Name the modified file: 'calc_fcc_ver1.in'. If pasting, use the 'Paste Special' command with 'Unformatted Text'. Note, this requires the variable latconst to be manually passed in from the command line when executing lammps, e.g. 'lmp_win_no-mpi -in calc_fcc_ver1.in - var latconst 4'.

# Find minimum energy fcc configuration
# Mark Tschopp, 2010
# This requires the variable latconst to be input via the command line
# e.g., lmp_win_no-mpi -var latconst 4 < calc_fcc_ver1.in

# ---------- Initialize Simulation --------------------- 
clear 
units metal 
dimension 3 
boundary p p p 
atom_style atomic 
atom_modify map array

# ---------- Create Atoms --------------------- 
lattice 	fcc ${latconst}
region	box block 0 1 0 1 0 1 units lattice
create_box	1 box

lattice	fcc ${latconst} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1  
create_atoms 1 box
replicate 1 1 1

# ---------- Define Interatomic Potential --------------------- 
pair_style eam/alloy 
pair_coeff * * Al99.eam.alloy Al
neighbor 2.0 bin 
neigh_modify delay 10 check yes 
 
# ---------- Define Settings --------------------- 
compute eng all pe/atom 
compute eatoms all reduce sum c_eng 

# ---------- Run Minimization --------------------- 
reset_timestep 0 
fix 1 all box/relax iso 0.0 vmax 0.001
thermo 10 
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 
min_style cg 
minimize 1e-25 1e-25 5000 10000 

variable natoms equal "count(all)" 
variable teng equal "c_eatoms"
variable length equal "lx"
variable ecoh equal "v_teng/v_natoms"

print "Total energy (eV) = ${teng};"
print "Number of atoms = ${natoms};"
print "Lattice constant (Angstoms) = ${length};"
print "Cohesive energy (eV) = ${ecoh};"

print "All done!"

Run this using LAMMPS in Windows

Follow these steps:
  1. Click on the Start button on the toolbar
  2. Click on Run...
  3. Enter 'cmd' and hit OK
  4. In new window, change to the directory that contains the LAMMPS executable (lmp_win_no-mpi.exe), the input script (calc_fcc_ver1.in), and the potential file (Al99.eam.alloy). For example, if these are on the desktop, type 'cd h:/desktop' to change to the desktop.
  5. Type 'lmp_win_no-mpi -var latconst 4 < calc_fcc_ver1.in'
  6. If the simulation completes, then it will type "All done"!

The end of the logfile/screen output should look like this:

Total energy (eV) = -13.43999995; 
Number of atoms = 4; 
Lattice constant (Angstoms) = 4.05; 
Cohesive energy (eV) = -3.359999988; 
All done!

Try modifying the starting lattice constant from the command line, e.g., to start with a lattice constant of 3.0 Angstroms, type 'lmp_win_no-mpi -var latconst 3.0 < calc_fcc_ver1.in'. This will not affect the final values (above), but it will start the simulation with a lower lattice constant for the FCC structure.

Run this from MATLAB

In this section, I will show how MATLAB can be used to run a LAMMPS simulation in Windows. Follow these steps:
  1. Open MATLAB
  2. Change to the directory that contains the LAMMPS executable (lmp_win_no-mpi.exe), the input script (calc_fcc_ver1.in), and the potential file (Al99.eam.alloy).
  3. Type 'system('lmp_win_no-mpi -var latconst 4 < calc_fcc_ver1.in');'
  4. The simulation should show up on the MATLAB screen. When it completes, it will display "All done"!

Now that you can run a simulation from MATLAB, it is easy to loop over a number of different lattice constants to calculate the energy. For instance, modify the input script again to only display the energy for the lattice constant input through the command line. The script should look like this:

# Find minimum energy fcc configuration
# Mark Tschopp, 2010
# This requires the variable latconst to be input via the command line
# e.g., lmp_win_no-mpi -var latconst 4 < calc_fcc_ver1.in

# ---------- Initialize Simulation --------------------- 
clear 
units metal 
dimension 3 
boundary p p p 
atom_style atomic 
atom_modify map array

# ---------- Create Atoms --------------------- 
lattice 	fcc ${latconst}
region	box block 0 1 0 1 0 1 units lattice
create_box	1 box

lattice	fcc ${latconst} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1  
create_atoms 1 box
replicate 1 1 1

# ---------- Define Interatomic Potential --------------------- 
pair_style eam/alloy 
pair_coeff * * Al99.eam.alloy Al
neighbor 2.0 bin 
neigh_modify delay 10 check yes 
 
# ---------- Run 0 --------------------- 
run 0

variable natoms equal "count(all)" 
variable teng equal "pe"
variable length equal "lx"
variable ecoh equal "v_teng/v_natoms"

print "Total energy (eV) = ${teng};"
print "Number of atoms = ${natoms};"
print "Lattice constant (Angstoms) = ${length};"
print "Cohesive energy (eV) = ${ecoh};"
print "%% ecoh = ${ecoh};"

print "All done!"

The 'run 0' command can be used to calculate the energy of the system without actually running an iteration. Name this 'calc_fcc_ver2.in' and combine this with a simple MATLAB script like this:

% MATLAB Script for running LAMMPS multiple times

count = 0;
for i = 3.0:0.10:5.0
   command_line = ['lmp_win_no-mpi -var latconst ' num2str(i) ' < calc_fcc_ver2.in'];

   % this next line executes the command line
   system(command_line) 

   % all that is left is to mine the 'log.lammps' file for the energy
   fid = fopen('log.lammps');
        tline = fgetl(fid);
        while ~feof(fid)
           matches = strfind(tline, '%%');
           num = length(matches);
           if num > 0 && matches == 1
                teval = strrep(tline,'%%','');
                eval(teval)
           end
           tline = fgetl(fid);
        end
    fclose(fid);

   % store the values in a matrix
   count = count + 1;
   X(count) = i; Y(count) = ecoh;
end

plot(X,Y,'-*r'), axis square

This script increases the lattice constant from 3.0 to 5.0 Angstroms in increments of 0.10 Angstroms AND calculates the energy for each lattice constant. As an easy extension, MATLAB can now be used to run hundreds of simulations, extract energies, and plot the lattice constant-energy curve in a matter of a minute or so. Notice the line 'print "%% ecoh = ${ecoh};"' in the 'calc_fcc_ver2.in' script - the '%%' is used as a flag for when MATLAB mines the log.lammps file. This MATLAB script iteratively reads every line in the log file and checks to see if the flag %% is contained within each line. When MATLAB finds a line with the flag, it removes the flag, and evaluates the rest of the line in MATLAB (which converts the string to a variable in the MATLAB stack).

Questions

Q: So, what does the 'run 0' command actually do? Is it minimization? Is it molecular dynamics?

A: The 'run 0' command is not minimization or molecular dynamics. The atoms are not displaced at all by this command. This is merely a way for LAMMPS to know that it needs to calculate additionally information for each atom or for the entire system. For instance, before the 'run 0' command, LAMMPS may not know the forces on each atom or the energy for the entire system. By executing the 'run 0' command, LAMMPS will calculate all the forces and attributes necessary to run molecular dynamics, but will not actually perform the timestep at the end, moving the atoms. Now, the user can print these attributes to screen or to a dumpfile, etc.

In some cases, the user will need to tell LAMMPS what values they want LAMMPS to calculate prior to a 'run 0' command. For example, if you want the stresses on every atom, you will need to define a compute or a variable command prior to 'run 0' and then let LAMMPS know that you need this value by inserting it in to the 'thermo_style' command (like in Tutorial 3!). For more information, see the LAMMPS documentation on the run command.