Problem with simulation of NaCl crystal

I removed the hybrid and fixed the neighbor skin. I am still getting the issue that is making my runs incredibly inconsistent with literature. In short, what happens is that a great amount of kinetic energy is added to my (supposed to be near fixed) atoms and destroys the lattice. In the LAMMPS documentation, the pair style equation appears with two changes from most literature:

  1. The variable 1/rho is used as B by many sources, so I just calculated B and took the inverse for the rho parameter.
  2. The D term appears as negative, where it appears as positive in the LAMMPS documentation. I just took the negative of the values given in the literature.

The issue is replicating itself.


Can you point out the source for that particular parametrization?
Would like to compare it to some old files I have for NaCl.


Sure thing I can. Recalling that I am in metal (and as such I will need eV for energy and A for length):

SOURCES: - Anwar et al. Calculation of the melting point of NaCl by molecular simulation. - Koishi et al. Large scale molecular dynamics simulation of aqueous NaCl solutions.

With that in mind, these are the parameters that they supply and my calculated ones:

Fumi-Tosi Units Aij (kJ/mol) B(1/A) Cij (A^6 * kJ/mol) Dij (A^8 * kJ/mol) Sigma_ij (A)

25.4435 3.1546 101.1719 48.1771 2.34

20.3548 3.1546 674.4793 837.077 2.755

15.2661 3.1546 6985.6786 14031.5785 3.17
LAMMPS Parameter Aij (eV) Rho (A) Sigma_ij (A) Cij (A^6 * eV) Dij (A^8 * eV)

0.2637055573 0.3169974006 2.34 1.0485818488 -0.4993247392

0.2109644458 0.3169974006 2.755 6.9905453129 -8.6757661783

0.1582233344 0.3169974006 3.17 72.4020778615 -145.4283108698

I’ve used that 1 eV = 96.4845 kJ/mol. Also note that as I mentioned, I had to invert the value of B (to obtain rho) value of D was made negative since their parametrization is positive, and LAMMPS is negative. I have tried both +D and -D, and it doesn’t allow me to converge to a stable rock salt crystal lattice at 298 K.

Thanks for any help you may be able to give regarding my parametrization,

To make it a little more clear, as I never placed who is what:

I think I found your problem. If I measure the distance between the Na and Cl ions in the lattice I am getting 5.something Angs. This feels way too large almost as if the units you used to set the lattice were
Bohr. Because your short-range cutoff = 7.6 < second-neighbor distance = sqrt(2)*5.something ~ 7.9
you artificially ended up with a cubic lattice that only has nearest-neighbor interactions. A cubic lattice with only NN interactions has zero shear resistance, i.e. is not stable, in other words does not correspond to a minimum of the potential energy surface. Just like shearing a cardboard box without the top and the bottom :wink:


Based on the literature, yes, you are right on the lattice parameter. I am using this equilibrium lattice distance to create my lattice actually. Also, your explanation makes complete sense to me and seems like the issue was pinpointed well by you.

I completely forgot about the NN distance, which would completely explain why I am having issues with my simulation. It’s funny that the literature I found says they are using this LJ cutoff distance. I will verify this data with my runs and get back to you replying whether everything is fine or not.

Thanks a million for your help.