Problem with simulation of NaCl crystal

Hello everyone,

Currently, I’m trying to simulate the melting of a sodium chloride (NaCl) crystal using the Born-Meyers-Huggins interaction with an added coulombic interaction term (via born/coul/long). I am generating the crystal within LAMMPS utilizing the commands:

/// ATOM GEOMETRY START ///

Na atoms

lattice fcc 11.2804 origin 0 0 0 # Lattice for Na system
region box block 0 2 0 2 0 2
create_box 2 box
create_atoms 1 box

Cl atoms

lattice fcc 11.2804 origin 0.5 0 0 # Lattice for Cl system
create_atoms 2 box

Other data

mass 1 28.990
mass 2 35.453
set type 1 charge +1.0
set type 2 charge -1.0

Misc.

neighbor 0.3 bin
neigh_modify delay 0

/// ATOM GEOMETRY END ///

And the generation is correct and no atom duplication occurs on the ends of the box (as checked per visualization in VMD). However, when I try to add velocity and time integration components to my system using a NVT integrator and velocity rescaling with the commands:

velocity all create 298.0 43454 dist gaussian mom yes
fix int all nvt temp 298.0 298.0 0.005

Quite inmediately, my NaCl crystal becomes highly disordered and while it does not melt, it actually does not follow a NaCl lattice at all either. Do you have any idea what could be ocurring that could be causing the issue? I am inclined to believe the problem is with the implementation of the parameters for the Born interaction term. However, these are taken (and compared) against a lot of sources in literature and seem to have no problem. Any help would be helpful and greatly appreciated. I am also attaching my implementation of the NaCl Born type interaction terms:

/// FORCE FIELD START ///

pair_style born/coul/long 7.6 15
pair_modify tail yes
pair_coeff 1 1 0.2637055573 0.317 2.34 1.0485818488 -0.4993247392 # Na-Na parameters
pair_coeff 2 2 0.1582233344 0.317 3.17 72.4020778615 -145.4283108698 # Cl-Cl parameters
pair_coeff 1 2 0.2109644458 0.317 2.755 6.9905453129 -8.6757661783 # Na-Cl parameters

kspace_style ewald 1.0E-5

/// FORCE FIELD END ///

Thanks,
Rafael Soler

I forgot to add, my idea is to first test the parameters agree with those at 298 K before melting (to ensure nothing out of order is happening).

Thanks,

Rafael,
What units are you intending to use? I don't see any units command in
your script. Note that the default units in Lammps are LJ which I'm
sure are not the ones you are trying to use.
Carlos

Hello everyone,

Currently, I'm trying to simulate the melting of a sodium chloride (NaCl)
crystal using the Born-Meyers-Huggins interaction with an added coulombic
interaction term (via born/coul/long). I am generating the crystal within
LAMMPS utilizing the commands:

# /// ATOM GEOMETRY START ///
        # Na atoms
lattice fcc 11.2804 origin 0 0 0 # Lattice for Na system
region box block 0 2 0 2 0 2
create_box 2 box
create_atoms 1 box
        # Cl atoms
lattice fcc 11.2804 origin 0.5 0 0 # Lattice for Cl system
create_atoms 2 box
        # Other data
mass 1 28.990
mass 2 35.453
set type 1 charge +1.0
set type 2 charge -1.0
        # Misc.
neighbor 0.3 bin
neigh_modify delay 0
# /// ATOM GEOMETRY END ///

And the generation is correct and no atom duplication occurs on the ends
of the box (as checked per visualization in VMD). However, when I try to
add velocity and time integration components to my system using a NVT
integrator and velocity rescaling with the commands:

velocity all create 298.0 43454 dist gaussian mom yes
fix int all nvt temp 298.0 298.0 0.005

Quite inmediately, my NaCl crystal becomes highly disordered and while it
does not melt, it actually does not follow a NaCl lattice at all either. Do
you have any idea what could be ocurring that could be causing the issue? I
am inclined to believe the problem is with the implementation of the
parameters for the Born interaction term. However, these are taken (and
compared) against a lot of sources in literature and seem to have no
problem. Any help would be helpful and greatly appreciated. I am also
attaching my implementation of the NaCl Born type interaction terms:

what units are you using? i don't see a units keyword anywhere, but your
parameters seem to hint at something else than the default of reduced (=
lj) units. also, the neighbor skin of 0.3 is inconsistent with the cutoffs.

axel.

Hi again,

Sorry for that (very huge) mess up. I’m trying to use metal units. The top part of my script (where I define the units) is:

/// DEFINITIONS START ///

dimension 3
units metal
boundary p p p
atom_style hybrid atomic charge

/// DEFINITIONS END ///

I’m verifying my units for the interactions to ensure I did not do any mistakes meanwhile.

Thanks,
Rafael

Hi again,

Sorry for that (very huge) mess up. I'm trying to use metal units. The top
part of my script (where I define the units) is:

# /// DEFINITIONS START ///
dimension 3
units metal
boundary p p p
atom_style hybrid atomic charge
# /// DEFINITIONS END ///

I'm verifying my units for the interactions to ensure I did not do any
mistakes meanwhile.

well, there is no need to use a hybrid atom style. the atom style charge
includes atomic.

axel.

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.

Thanks,