Questions about Ni-Si alloy with the MEAM potential

Hi all,

I have been having some issues with the MEAM potential in LAMMPS. I have been attempting to use MEAM to model Ni-Si interactions, and have been using the implementation discussed by Baskes et al in the paper “Atomistic calculations of composite interfaces” in the journal Modelling Simul. Mater. Sci. Eng. 2 (1994) 505-518.

According to Greg Wagner in a previous post I made (, the “l12” structure is properly coded in LAMMPS and the Ni-Si interactions can be easily modeled by using the parameters given in the paper. Previously, I have been using their parameters to model amorphous-Si/Ni interactions, the structure of which deviates significantly from the reference structures used to define the interactions. Although the behavior of the Ni-Si interactions has been questionable, I thought that that would be primarily due to the use of the Baskes Ni-Si parameters for a system far from the reference structure. However, it appears that if I use the Baskes parameters, I cannot even correctly describe the reference structure from which they were defined (Ni3Si). Can anyone tell me if I am doing something wrong or why the Baskes Ni-Si parameters are unable to model the reference Ni-Si structure?

I’ll include a simple input script that is just a block of Ni3Si sitting at 300K. The parameters I’m using match those in the Baskes paper, but the system quickly loses its crystal structure (regardless of the equilibration) and goes unstable (rapid rise in temperature and energy) depending on the length of equilibration. But a pure Ni (fcc) and pure Si (diamond) block maintain their structure and stability without any equilibration. Any ideas would be greatly appreciated.

Michael Price

P.S. I have been using the 9-Nov-13 version of LAMMPS (behavior is the same using earlier versions as well)

LAMMPS input script:


variable outputfolder string output
shell mkdir ${outputfolder}

dimension 3
boundary p p p
units metal
atom_style atomic
neighbor 2.0 bin
neigh_modify delay 5
lattice fcc 3.405

region box block 0 10 0 10 0 10
create_box 2 box

mass 1 58.71 # Ni
mass 2 28.086 # Si

pair_style meam
pair_coeff * * library.meam Ni1 Si2 NiSi.meam Ni1 Si2

#-------------CREATE ATOMS-----------------------------------

create_atoms 1 region box basis 1 2 basis 2 1 basis 3 1 basis 4 1


velocity all create 300 1272
fix int all nve
timestep 0.0001

thermo 50
thermo_style custom step temp epair emol etotal press
dump 1 all custom 100 ${outputfolder}/dump.stats id type x y z
dump_modify 1 sort id

dump 2 all cfg 100 ${outputfolder}/conv*.cfg mass type xs ys zs
dump_modify 2 sort id
dump_modify 2 element Ni Si


System quickly goes unstable unless there is a fairly long (>10ps)

equilibration process and small (.00025 ps) timestep. A pure

Ni (fcc) or pure Si (dia) lattice is stable even without


fix equilibrate all temp/rescale 1 300 300 1 1
#fix equilibrate all langevin 300 300 .01 39857

run 12000

unfix equilibrate

run 1000

--------------------------------------------MEAM FILES----------------------------

meam data from vax files fcc,bcc,dia 11/4/92

elt lat z ielement atwt

alpha b0 b1 b2 b3 alat esub asub

t0 t1 t2 t3 rozero ibar

‘Ni1’ ‘fcc’ 12. 28 58.71
4.99 2.45 2.20 6 2.20 3.52 4.45 1.10
1.0 3.57 1.60 3.70 1.0 0
‘Si2’ ‘dia’ 4. 14 28.086
4.87 4.4 5.5 5.5 5.5 5.431 4.63 1.
1.0 3.13 4.47 -1.80 1 0


rho0(1) = 4.88
rho0(2) = 1
delta(1,2) = 0.36
alpha(1,2) = 6.00
re(1,2) = 2.41
lattce(1,2) = l12

I’d say this is a Q for Greg (CCd).


Hi Michael,

As a first step, have you verified that the initial configuration (which should just be the unperturbed L12 lattice) gives the correct energy per atom? If not, then there’s something wrong with either the parameters or the atom placement.

I haven’t verified your input file, but a few things I’d check:
— Are the atoms really being placed with correct initial spacing?
— It seems strange that you don’t have a cohesive energy (Ec(1,2)) defined in the parameter file.
— I’m not sure about this, but you may need to have the lattice style ‘l12’ in quotes in the parameter file (it is in all the local examples I’ve looked at, and it’s been a while since I’ve played with this).

Attached are some files I used a few years ago to test a similar L12 system. The current run just sets up the configuration and computes the total energy. Note that the initial energy of 1090.688, divided by 256 atoms, matches the cohesive energy per atom of 4.2605 defined in the input file. Hopefully by comparing this input with your own system you’ll find something amiss.


Greg Wagner
Manager, Thermal/Fluid Science and Engineering Department
Sandia National Laboratories, Livermore, CA
Tel: (925) 294-2180 Fax: (925) 294-3410
Email: gjwagne@…3…

Fe3Cu_in.meam (369 Bytes)

Fe3Cu_param.meam (426 Bytes)

Fe3Cu_position.meam (7.3 KB)

library2.meam (9.54 KB)

log.lammps (1.17 KB)

Hi Greg,

Thanks for your suggestions and the Fe3Cu files. I have gotten the Ni3Si system working, although it has brought up a couple of questions.

First, in response to your suggestions: the atoms were being placed with the correct initial spacing, the single quotes for the ‘l12’ structure are actually not necessary, and the cohesive energy is defined as “(3E_Ni+E_Si)/4 - delta” if delta (heat of formation) is given but Ec is not. That is how it was given in the Baskes paper I am referencing, although I have not compared this calculated value to cohesive energies given in the literature. The (total energy) / (number of atoms) printed at timestep 0 does match the calculated cohesive energy, equivalent to what I saw with your code.

There were three things I did to get the Ni3Si system working, but if I’m correct then I may have found an error in the code. According to the LAMMPS documentation for the mair_style meam, “The ibar parameter selects the form of the function G(Gamma) used to compute the electron density”, and various options are given. Looking at the code (meam_setup_done.F line 353 in the most recent version), the equations coded in LAMMPS reference “I. Huang et al, Modelling Simul. Mater. Sci. Eng. 3:615”, which uses the ibar = 1 form of the G_gam function. In that paper, the equation is exp(-gamma/2) but in the code it is implemented as exp(gamma/2), without the negative sign (meam_dens_final.F lines 222 and 267). Is that a mistake or is that intentional? I cannot get the Ni3Si system to work using the code as written, but if I include the negative sign in the source code, then it works. The inclusiong/exclusion of the negative sign has negligible effect on simulations of pure Ni, pure Si, or the Fe3Cu code you sent me, but does affect the Ni3Si.

The the Ni3Si system works with ibar = 1 (with the negative exponential) but does not work with ibar = 0. This seems strange considering that the Baskes parameters I am using come from a paper using the ibar = 0 form. Do you know if that’s reasonable or if it’s an indication that I am doing something wrong?

Finally, I had to change the rho0 value for Ni from 4.88 to 1 and include the Cmin and Cmax terms from your Fe3Cu code. I saw that in your Fe3Cu code there was no relative density scaling, so I tried it with Ni3Si and it worked. It appears that the implementation in the code takes care of the relative scaling, so the value given in the Baskes paper (4.88) is not necessary (and actually makes it unstable). The Cmin and Cmax parameters are also necessary to get Ni3Si to work, but the values from your Fe3Cu code are not optimal for Ni3Si - without the Cmin and Cmax values, the system is unstable, but with the Cmin and Cmax values the (total energy) / (number of atoms) printed at timestep 0 now has an error of about 1% from what it should be. Do you have any suggestions or references I can use to get better Cmin and Cmax values for Ni3Si?

Again, the Ni3Si system seems to be working with the changes I’ve mentioned. Thanks again for your initial comments and help, they have helped a lot. I would greatly appreciate any other advice or comments you have.


Hi Michael,

Thanks for looking so deeply into the MEAM code. I’m sorry it’s been giving you trouble. For reasons I’ll describe I don’t think there’s a sign error in the G_gam function, but unfortunately I don’t have too many more insights that might be helpful right now.

The options for the G_gam function come from Baskes’ paper, “Determination of modified embedded atom method parameters for nickel,” Materials Chemistry and Physics 50 (1997) 152-158. Equation (12) has the expressions used for G(Gamma), and (12b), which is the relevant one here, has a positive sign for Gamma. Maybe more importantly, note that all of the forms of G(Gamma) are increasing for increasing values of Gamma, which wouldn’t be true if Gamma were negative in the expression. In fact, as this paper points out all of the forms of G have the property that G’(0) = 1/2, which again requires a positive exponent here.

I can’t explain the Huang paper — I assume you’re looking at equations (6) and (7). It could be a typo, but it’s annoying that it’s repeated twice. The only other data point I’ll add is the paper by Baskes and Johnson from the previous year (MSMSE 2: 147-163), which in equation (4c) has the exponential form without the minus sign.

I realize that none of this is helpful if it’s not giving you stable results. I don’t see anything else wrong with your inputs. Have you verified any of the computed quantities from the MSMSE paper you’re taking your parameters from? If those match, then it’s very possible that the potential in the paper truly does give unstable results. In the paper, I only see static quantities compared, and much of the focus is on the Ni-Si interface (not the L12 structure).

You might want to e-mail Mike Baskes directly about this. He might have some insights, and may know whether this potential was ever tested on dynamic simulations.


Greg Wagner
Manager, Thermal/Fluid Science and Engineering Department
Sandia National Laboratories, Livermore, CA
Tel: (925) 294-2180 Fax: (925) 294-3410
Email: gjwagne@…3…