Incorrect well-depth for lj/cut potential

Dear LAMMPS mailing list,

In the framework of my project, I am using the Lennard-Jones potential to describe the argon – silicon interactions for the simulation of argon bombardment of silicon. By fitting the Lennard-Jones potential 12280404.gifto DFT data (Figure joined in annex), I have obtained the epsilon and sigma parameters for the argon – silicon interactions. The parameters give a well depth of 0.02298 eV (0.52854 kcal/mol). I am using the 5 May 2020 GitHub patch but this problem also appeared in previous LAMMPS versions.

In order to check the correctness of my parameters in LAMMPS, I output the potential energy as a function of distance for a pair of Ar-Si atoms (script in annex). To output the energy of my atoms, I am using the compute pe/atom and compute ke/atom commands. I run this script and change the distance between particles by increments, outputting the results of each compute in a file in order to reconstruct the data shown in the figure joined in annex, representing the obtained potential between my particles. When doing so, I can only tune two parameters of the lj/cut pair_style, which are the well depth and the zero crossing distance. After putting the aforementioned parameters into LAMMPS I obtain a shifted zero crossing distance and my depth well is divided by a factor 2!

In fact, to obtain the same well depth than the DFT data, I have to put 1.05712 kcal/mol for the well depth. Even tho the zero crossing distance seems to be shifted, the position of the potential’s minimum is correct, and it is the most relevant information for my case studies therefore I am not worrying too much about it.

My question is therefore: how are the well depth and the zero crossing distance handled by the lj/cut command, and why do I have to multiply by two the well depth in order to reproduce my DFT data? In the case the DTF would be wrong, the reproduced values in LAMMPS should still be fitting the DFT data, yet I have this factor two that I cannot explain on my own or after discussing it with my Ph.D. supervisor. Maybe the lj/cut command is not suitable for my specific case ? Anyone has already faced this issue ? Any help would be greatly appreciated.

Thank you in advance,
Best regards,
Grégoire DEFOORT, Ph.D.

Figure of the potential outputs:

(See attached file: Figure_Potentials.png)

Script Used:

3d reaxFF system

— initialisation ----------------------------

units real
#Real units: Kcal/mole - Angstroms - femtoseconds
dimension 3
boundary f f f

atom_style full
#atom_modify map array

print “Initialisation done…”
print “”

— variables ---------------------------------

variable time_step equal 1e-4
variable m_si equal 28.0855
variable m_ar equal 39.948

— atoms definition --------------------------

lattice none 1
region box block -10 10 -10 10 -10 10 units box
create_box 2 box

region silicon block -10 10 -10 10 -10 10 units box

#Fixed position Silicon particle
create_atoms 1 single 0 0 0

#Shifted position Argon particle
#The script loops and increases the z value of the Ar particle by 0.01 each step
create_atoms 2 single 0.0 0.0 0.01

print “Atoms defined”
print “”

— settings ----------------------------------

pair_style lj/cut 10
pair_coeff 1 2 1.05712 3.624 10
pair_coeff 1 1 0.98899999609 3.428 10
pair_coeff 2 2 0.69046 3.40776 10

mass 1 {m_si} mass 2 {m_ar}

neighbor 2.5 bin
neigh_modify every 1 delay 0 check yes

timestep ${time_step}
run_style verlet

— computes & thermo -------------------------

compute pe_Atom all pe/atom
compute ke_Atom all ke/atom

compute total_potential_energy all pe
compute total_kinetic_energy all ke

thermo 100
thermo_style custom step dt time temp density pe ke etotal press

dump 1 all custom 1 energy_Si-Ar.nve id element c_ke_Atom c_pe_Atom
dump_modify 1 element Si Ar format float “%20.13f”

thermo_modify lost warn norm yes flush yes

print “Settings done”

— run ---------------------------------------

print “Starting simulation”
run 1
print “Simulation ended”



i don’t want to discuss here the validity (or problems) of DFT for interactions which include dispersion, and thus the validity of using DFT for your input data here. that is a separate topic and has nothing to do with LAMMPS and thus outside of the scope of this mailing list.

also, i think we can safely assume that the lj/cut potential in LAMMPS is implemented correctly. LAMMPS has been around for 25 years this year and somebody would have noticed, given how much it is used and that its results match those of other MD codes.

thus the cause of the unexpected outcome of your potential fitting has to be a result of using an inadequate fitting procedure. as you can easily see from your graphs the shape of your DFT data does not fit the LJ potential very well. especially the repulsive part of the 12-6 LJ is way too steep. this is a known property of the 12-6 LJ and that leads to problems when applying a least squares fitting procedure (which i assume you are using). you can verify this by using different fitting ranges. you should get significantly different parameters from the fit depending on the range of distances you use. thus a more meaningful fit should be obtained from fitting to the location of the minimum (which is at sigma*2^(1/6)) and the curvature near the minimum.

I would suggest to also redo the fit using a morse potential, which has a softer and more physically correct repulsion. the 12-6 LJ exponents are chosen for computational efficiency not physical correctness, but the impact of that is usually limited since most of the time the interactions are around the minimum or in the attractive branch, where the match is not so bad.
after obtaining a (likely much better) least squares(?) fit to the morse potential you could then infer LJ parameters analytically from the fitted morse potential by matching well depth and minimum. that would be a more consistent approach. however, you cannot get around the mismatch of the shape between your DFT data and a 12-6 LJ, so you will have to make a choice on which range of the potential function you want to match the DFT data best.
or just use the morse pair style right away. there are other pair styles available in LAMMPS that may give an even better representation.

please let me also point out that LAMMPS has a “table” pair style, which can represent any shape of a potential, you would only have to obtain the forces from numerical differentiation (which can be done rather crudely through finite differences or through more sophisticated methods as from a spline fit or some piece-wise fitting of a higher order polynomial and then taking the derivative of that polynomial, aka as Savitzky-Golay filtering )




Hello Grégoire,

Also note that LAMMPS has a builtin pair_write function which computes the force and energy between pairs of atoms and uses the tabulated format as an output. Overenginneering using such a complicated procedure is error prone.
Moreover, the attached script does not do what you say it does (no position update) so can’t tell what you actually plotted. With a quick test using the pair_write command I got the curve you expected from your fitted parameters with correct depth and minimum.

As Axel mentionned lj/cut is a commonly used and well tested potential in LAMMPS so it is very unlikely there are any error in the computation of energy or forces. Check with your advisors what you are actually computing with regards to pair interaction energy.


Dear LAMMPS mailing list,




for Argon sputtering of Si, using ZBL or Molière potential would Be a good approximation. Your LJ potential is giving a harder repulsive part than DFT while it is the most important one. 0.022 eV is below room temperature and any high velocity Ar impinging will not feel the depth well. Moreover the depth well will be reduced to thermal transfer if Si is at 300K and above. Depth will be meaningful at very low Ar kinetic energy and below Si surface temperature 190 K, temperature at which Ar could bind on a surface. Summarily, it makes almost no sense to search for Ar-Si,not well fitting LJ parameters for such a problem.
I would not recommend such approach.

Pascal Brault
CNRS-Université d’Orléans

Dear Pascal, Germain and Axel,

First of all, thank you for your very interesting replies. Maybe my words were misplaced / not clear enough, as I know how well defined the lj/cut function is and I didn’t want to say that it is badly coded, but more that the error was within my end of the coding.

I will definitively dig the pair_style table and pair_write function in order to simplify my code and verify my output, it will help my situation. Indeed, after checking the litterature, a Morse potential could be usefull to estimate the LJ parameters in my case study. The savitzky - Golay filter is a vald methodology that I will apply to my DFT data. I am affraid tho that by using a table for my potential, some interpolation have to be made (in opposition to a well defined potential), and could cause some discrepancies over a lot of steps with some error margins. I am in fact studying extremely low energy sputtering, so the depth well (even tho overlapped by thermal transfers, as stated by Pascal) is crutial to my case study.

Thank you again for you inputs and comments,
Best regards


graycol.gif“Pascal Brault” —02/06/2020 13:28:43—Hi Moreover,

Dear Pascal, Germain and Axel,


I am affraid tho that by using a table for my potential, some interpolation have to be made (in opposition to a well defined potential), and could cause some discrepancies over a lot of steps with some error margins.

Sorry, but that argument makes not sense to me. If that would be a concern, then you must not use a classical MD method at all, but do your simulation with DFT.
A tabulation, if done properly, can adapt better to your DFT data, which is tabulated data to begin with, than a fit to a specific shape of a given functional form. Your graphs clearly show that there is a significant mismatch, specifically for the repulsive branch.