[lammps-users] Electrostatic interaction in NaCl

Hello everyone,

I am trying to simulate the structure of different inorganic structures (e.g. silica, rutile) via interatomic potentials (e.g. BSK potential). Whatever potential I use, I found unrealistic pressure values that are about -10**5 bar and a Coulomb potential two order of magnitude bigger than the kinetic energy and other VdW potentials.

I can reproduce this problem (unrealistic pressure and large negative value of the Coulomb potential) with a NVT simulation of crystal NaCl at 300K, so I will use this system as an example. The electrostatic energy (E_coul + E_long) computed for this system is -35.7 eV. As far as I understand, this value should correspond in magnitude to the Madelung potential of NaCl, which is roughly 9 eV. What is wrong with my input files? I am attaching them at the end of the email.

Thank you very much for the help,


Are you normalizing by the number of atoms?


Hi Steve,

Are you normalizing by the number of atoms?

Touché! As the (normalized) coulomb interaction for NaCl is correct, I move to the real problem.

I want to use the BSK potential [PRL 64, 1955 (1990)] to simulate the structure of amorphous silica. This potentail is written as the sum of the Coulomb and Buckingham potentials, and it has no Si-Si interaction.
I wrote it in the following way:

# Si type 2, O type 1
kspace_style pppm 1.0e-4
pair_style buck/coul/long 10. 12.
pair_coeff 1 1 1388.7730 0.362319 175.0
pair_coeff 1 2 18003.7572 0.205205 133.5381
pair_coeff 2 2 0.0 1.0 0.0

I have simulated a box (57.8 x 57.8 x 62.6 A**3) of SiO2 in the NVT ensemble at 300K using this potential and a supercell of cristobalite, and obtained a very large negative pressure (~ -24000 bar). For this system I expect the pressure to be around zero, as the starting structure has the correct density at RT. It seems to me that attractive forces are too strong, and as a consequence atoms are often lost during the trajectory. Is there something missing in the specification of the BSK potential?

While waiting, I am doing further tests on silica using the Tsuneyuki potential [PRL 61, 869 (1988)].

Thank you very much for the help, I do appreciate it.

I don't know. I would start with a simpler small bulk system and
see if you can reproduce thermodynamic state info in the literature.