LJ cutoff question

Dear LAMMPS users,

When I’m trying to do a simple test using the 12-6 type LJ potential, I found the results are very sensitive to the cutoff specified in the pair_style command. I know it is a convention to use a LJ cutoff distance of 2.5 times sigma, in my case sigma=3.145 A. However, even if I set the cutoff much larger (ranging from 10 A to 40 A), after energy minimization, the lattice constant and the minimum energy per atom still changes. I have attached my input script and listed my results below. Can anyone please help give me some suggestions or insights?

cutoff(A) lattice constant(A) energy/atom(eV)
10 6.325 -11.032
20 6.4889 -13.437
30 6.6667 -14.183
40 6.7123 -14.22

Thanks very much!

Chun Tang
University of California

in.lj (652 Bytes)

Dear LAMMPS users,

When I'm trying to do a simple test using the 12-6 type LJ potential, I
found the results are very sensitive to the cutoff specified in the
pair_style command. I know it is a convention to use a LJ cutoff distance of
2.5 times sigma, in my case sigma=3.145 A. However, even if I set the cutoff

that is primarily applicable for MD and 2.5 is a quite aggressive
lower limit. 3.5-4.5 is more common these days. at finite temperature,
fixed volume and for a homogeneous system, the error from cutting off
the forces on the individual atom largely cancels.

much larger (ranging from 10 A to 40 A), after energy minimization, the
lattice constant and the minimum energy per atom still changes. I have
attached my input script and listed my results below. Can anyone please help
give me some suggestions or insights?

just think about this for a bit.
- the LJ potential is attractive for r > sigma**(1/6)
- the larger the cutoff, the more atoms you include (grows with O(r**3))

=> even for very large cutoffs you have non-negligible energy contributions.
to some degree it is comparable to long range electrostatics.

have a look at the documentation for the lj/long pair_style for use
with pppm/disp or ewald/disp

axel.

If you are doing energy minimization, are you using
the pair_modify shift yes option? If you are
not, then there is a jump in your energy at the
LJ cutoff which means the minimization is ill-defined.
The minimize doc page warns you about this.

Steve

Dear Steve and Axel,

Thanks very much for your kind instructions, they are very helpful!

I’m now having a question related to this: it is known from textbooks that the cut off should be no more than half the size of a unit cell in simulation so as to avoid self interactions or interacting with the same atom multiple times. In a recent test, I however, found that using a cutoff of 50 A and cell sizes ranges from 6.29 to 125.8 A, the obtained elastic constants are the all same: 248 GPa. Are they supposed to be different?

The potential I used is a tabulated table in combination with coul/long potential for the KCl system. I’m attaching one of my input examples and potential here, could anyone please advise why this is so?

Thanks very much!

Chun

init.mod (1.49 KB)

in.elastic (5.79 KB)

potential.mod (821 Bytes)

displace.mod (4.51 KB)

BMH_50a.lammps (313 KB)

Dear Steve and Axel,

Thanks very much for your kind instructions, they are very helpful!

I'm now having a question related to this: it is known from textbooks that
the cut off should be no more than half the size of a unit cell in
simulation so as to avoid self interactions or interacting with the same
atom multiple times. In a recent test, I however, found that using a cutoff
of 50 A and cell sizes ranges from 6.29 to 125.8 A, the obtained elastic
constants are the all same: 248 GPa. Are they supposed to be different?

minimum image conventions do not apply to LAMMPS, since it creates
real copies (as "ghost atoms").

axel.

Dear Axel,

Thanks again for your kind help!

Regarding the kspace summation of long range energy. I found that whether or not using the following commands makes very big difference on the elastic constants of a crystal structure:

kspace_modify gewald 0.29

I have read the manual carefully but found how this number 0.29 is chosen is not instructed, could you please advise on how to estimate the value for gewald parameter?

Thanks very much!

Chun

Dear Axel,

Thanks again for your kind help!

Regarding the kspace summation of long range energy. I found that whether or not using the following commands makes very big difference on the elastic constants of a crystal structure:

kspace_modify gewald 0.29

I have read the manual carefully but found how this number 0.29 is chosen is not instructed, could you please advise on how to estimate the value for gewald parameter?

you have to read up on ewald summation. the LAMMPS manual explains how
to use LAMMPS, it is not meant to explain the underlying models and
methods. for that you have to refer to text books and research
articles. many of the latter are referenced at the appropriate places,
for example the article describing the method that LAMMPS uses to
estimate the kspace convergence.

axel.