We are hoping to use the new adiabatic core-shell model feature available in the 15May15 version of LAMMPS and found a descrepancy for a test system which may indicate that we are using the code incorrectly or may indicate that there is an energy shift in the core-shell code. We used LAMMPS to calculate the total potential energy of a periodic cubic unit cell of MgO using a buckingham potential with and without the core-shell feature. In this case, there should be no physical effect of the core-shell
because of the high symmetry, but LAMMPS calculates a small (~0.15 eV) shift in the potential energy. The results for the potential energy are
No core-shell “nocs” With core-shell “cs”
-159.5507 eV -159.40695
Further calculations at other lattice constants suggest that the shift may be constant and therefore perhaps inconsequential, but thought we should check with the experts. The input files for the two cases are attached. Thanks in advance for any advice. Natalie Holzwarth and Ahmad Al-Qawasmeh
in.lammps.nocs (560 Bytes)
readdata.nocs (948 Bytes)
in.lammps.cs (671 Bytes)
readdata.cs (1.77 KB)
Dear all,
The energy descrepancy is caused due to a non-optimal parametrization of the real-space part of the Ewald summation for short distances. Alain Dequidt who was involved in the USER-DRUDE package supplied a better parametrization. On basis of this we recreated the pair_cs styles including the improved parametrization as well as the already present internal minimal distance (enabling computation for core/shell distances of r=0 and eliminating oscillation of the real-space part of the Ewald summation for distances smaller and around r=1e-6 Anstroem). The new pair_styles have just been released in the new version of LAMMPS 9 Sep 2015.
You can see below that the energy deviation is now negligibly small (only caused by the internal minimal distance). If you also compare the forces (Fnorm and Fmax) you can see that the descripancy here is caused by very small (e-14) round-off errors proportional to the number of particles in your system (the force should normally be zero in your high symmetry case).
Step TotEng PotEng KinEng E_pair E_vdwl E_coul E_long E_bond Fnorm Fmax
no_cs 0 -159.5507 -159.5507 0 -159.5507 11.780834 -132.75149 -38.580044 0 3.2654059e-14 1.6209256e-14
cs_old 0 -159.40695 -159.40695 0 -159.40695 11.780834 -108.11692 -63.070861 0 8.6513875e-14 6.367129e-14
cs_new 0 -159.55064 -159.55064 0 -159.55064 11.780834 -108.26061 -63.070861 0 8.7918759e-14 5.1514348e-14
I also noticed that you created core/shell pairs for the non-polarized ions. Due to another recent mailing-list question (lammps.sandia.gov/threads/msg55255.html) I realized that we did not state explicitly in the documentation that you can mix polarized and unpolarized ions in the LAMMPS-Implementation of the core/shell model. At this point I would recommended that you remove the sattelite particles as this is computationally cheaper.
Best wishes