Use MEAM in MD++
Keonwook Kang, Seunghwa Ryu and Wei Cai
This tutorial describes how to use the modified embedded-atom method (MEAM) potential in MD++. The MEAM potential was originally developed by Dr. Mike Baskes,
Modified embedded-atom potentials for cubic materials and impurities, Physical Review B, 46, 2727 (1992).
Since then MEAM potentials for many elements, as well as binary and ternary alloy systems have been developed. The MEAM potentials allows us to quickly set up a model for a wide range of materials systems. MEAM has also been implemented in LAMMPS, whose subroutines have been used in the MEAM implementation in MD++.
Through some examples, we will describe how to use MEAM in MD++ to model a binary system, as well as allowing the simulation cell to change its shape.
Contents |
Download files
First, you need to install MD++ on your computer by following the instructions on MD++ Manuals.
Second, copy the input file meam-si-c.tcl and potential files meamf and SiC.meam to your input file directory of MD++. You can do so by the following commands (assuming you have installed MD++ in ~/Codes/MD++).
export MDPP=~/Codes/MD++
mkdir -p ${MDPP}/scripts/work/meam
cd ${MDPP}/scripts/work/meam
wget http://micro.stanford.edu/mediawiki-1.11.0/images/Meam-si-c.tcl.txt -O meam-si-c.tcl
wget http://micro.stanford.edu/mediawiki-1.11.0/images/Meamf.txt -O meamf
wget http://micro.stanford.edu/mediawiki-1.11.0/images/SiC.meam.txt -O SiC.meam
Compile executable file
There are several MEAM implementations in MD++,
meam-lammps - MEAM (taken from lammps) meam-baskes - MEAM (taken from Baskes's code dynamo) meam - MEAM (directly implemented in MD++)
meam-lammps is the one that is most extensively tested. We will use meam-lammps exclusively in the following discussions.
You can compile meam-lammps by typing
cd ${MDPP}
make meam-lammps build=R SYS=su-ahpcrc
If compilation is successful, this will create binary file meam-lammps_su-ahpcrc file in your ${MDPP}/bin directory. (su-ahpcrc is one of our computer clusters at Stanford.)
meam-lammps requires Fortran compiler. This is because it uses a number of Fortran codes taken from LAMMPS (in ${MDPP}/Fortran/MEAM-Lammps). Make sure that in the ${MDPP}/src/Makefile.base file, you have specified the correct Fortran compiler on your computer (e.g. gfortran on Mac), as well as the path of Fortran library (e.g. -L/usr/local/lib -lgfortran on Mac).
MD of Si crystal in NVT ensemble
The following command will run NVT ensemble for a perfect crystal made of 1000 Si atoms.
cd ${MDPP}
bin/meam-lammps_su-ahpcrc scripts/work/meam/meam-si-c.tcl 0
This is assuming you have put the si-meam.tcl input file in your ${MDPP}/scripts/work/meam directory. An X-window should open and display the positions of the Si atoms during the MD simulation. Output of the simulation is saved in the prop.out file. The variables saved in this file is specified by the output_fmt variable in the input file. For more details about finite temperature simulations in MD++, please read manual 05.
MD of Si crystal in NσT ensemble
The following command will run NσT ensemble for a perfect crystal made of 1000 Si atoms.
cd ${MDPP}
bin/meam-lammps_su-ahpcrc scripts/work/meam/meam-si-c.tcl 1
A shear stress σxy = 1000 MPa is applied to the simulation cell. All other stress components are zero. The simulation cell is allowed to adjust its shape during the simulation. Three of the nine component of the H matrix (whose three column vectors are the repeat vectors) are fixed (specified by fixboxvec) to remove arbitrary rotation.
The 7,8,9 columns of the output file prop.out saves the instantaneous values of H12, H22, H32 (in angstrom). We can see that H12 changes from 0 to -0.43 in response to the applied shear stress. (The negative sign is due to the sign convention of the stress in the Parrinello-Rahman boundary condition, J. Appl. Phys. 52, 7182, 1981.) This corresponds to a shear modulus C44 of approximately 60 GPa (at T = 1000 K). We may compare this with experimental data (80 GPa at T = 300 K). H22 changes from 27.2 to 27.5 due to thermal expansion.
The 10,11,12 columns of the output file prop.out saves the instantaneous values of σxy, σyy, σzy (in MPa). We can see that σxy eventually fluctuates around 1000 MPa as it should. The other stress components fluctuate around zero after equilibrium is reached.
After the MD simulation, the structure is relaxed by conjugate gradient algorithm with the box fixed.
MD of SiC crystal in NσT ensemble
The following command will run NσT ensemble for a perfect SiC crystal made of 500 Si and 500 C atoms.
cd ${MDPP}
bin/meam-lammps_su-ahpcrc scripts/work/meam/meam-si-c.tcl 2
This simulation requires not only the standard meamf file, but also the SiC.meam file to specify the interaction between Si (species = 0) and C (species = 1) atoms. augt1 = 1 in SiC.meam (for LAMMPS code) is equivalent to legend = 0 in meafile (for Baskes's dynamo code). This corresponds to the earliest version of MEAM potential where regular polynomials (instead of Legendre polynomials) are used. Other simulation settings are the same as the previous test case. From this simulation, we can estimate the shear modulus C44 (at T = 1000 K) to be approximately 200 GPa. We may compare this with ab initio result (237 GPa at T = 0 K).
The SiC/MEAM potential used in this case study is described in
Huang, Ghoniem, Wong and Baskes, Molecular dynamics determination of defect energetics in β-SiC using three representative empirical potentials, Model. Sim. Mater. Sci. Eng. 3, 615 (1995).
Technical Note
Because MEAM potential is quite complex, different implementations (e.g. MD++/LAMMPS and DYNAMO) can have small differences.
For example, in DYNAMO, the equation of state function contains the following line,
erose=-esub*(1.+x+an3*x**3/(r/re))*exp(-x)
where an3 = attrac (if x > 0) and an3 = repuls (if x < 0). x = alpha*(r/re - 1.0).
In MD++/Fortran/MEAM-Lammps/meam_setup_done.F, we can find the following line,
c erose = -Ec * (1 + astar + a3*(astar**3)/(r/re)) * exp(-astar)
erose = -Ec*(1+astar+(-attrac+repuls/r)*(astar**3))*exp(-astar)
where astar is the same as x in DYNAMO.
The first line should be the same as that in DYNAMO, but it is commented out. The second one has a different form for attrac (not divided by r). The two lines are equivalent only when attrac = 0, and when repuls is scaled by re(nearest neighbor distance). Of course, you can always modify the source code to make the equation of state function behave the way you want.
DYNAMO also has some special treatment when square of the electron density 
becomes negative. This unphysical situation occurs when two atoms become too close to each other. MD++/LAMMPS does not have this special treatment so that the numerical results can be different from DYNAMO.