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:
 The variable 1/rho is used as B by many sources, so I just calculated B and took the inverse for the rho parameter.
 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,
Rafael,
Can you point out the source for that particular parametrization?
Would like to compare it to some old files I have for NaCl.
Carlos
Carlos,
Sure thing I can. Recalling that I am in metal (and as such I will need eV for energy and A for length):
SOURCES:
http://jcp.aip.org/resource/1/jcpsa6/v118/i2/p728_s1  Anwar et al. Calculation of the melting point of NaCl by molecular simulation.
http://www.iapws.jp/Proceedings/Symposium03/212Koishi.pdf  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:
FumiTosi 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,
Rafael
To make it a little more clear, as I never placed who is what:
Rafael,
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 shortrange cutoff = 7.6 < secondneighbor distance = sqrt(2)*5.something ~ 7.9
you artificially ended up with a cubic lattice that only has nearestneighbor 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
Carlos
Carlos,
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.
Rafael