Truncated vs Truncated-Shifted LJ

Dear Lammps Users,

This is a question that has been discussed many times before – but I could not find the answer that I am specifically looking for.

What LAMMPS says:
since derivatives are identical (if I cut off anywhere above sigma), the trajectories and forces are identical. Relevant details from LAMMPS manual on pair_modify shift yes – “If so, this adds an energy term to each pairwise interaction which will be included in the thermodynamic output, but does not affect pair forces or atom trajectories.”

My Simulations:
I performed two simulations (one with truncated and other with truncated/shifted LJ, i.e., with and without pair_modify shift yes at 2^(1/6) sigma). I understand there is discontinuity in potential and force for case 1 and latter for case 2. Now to the question. I printed intensive evdwl from thermodynamic output with both the cases, and it gives identical results. Why is it so? Shouldn’t the extra offset be reflected in this or am I missing something. Relevant code details and outputs are give below. Please let me know if you need the full code and I shall provide the same.

Pair Information

bond_style fene
special_bonds fene
bond_coeff 1 30 1.6 1.5 1

angle_style cosine
angle_coeff 1 3
dihedral_style none
improper_style none

pair_style lj/cut/coul/long 1.122462 10.0
pair_coeff * * 1 1
kspace_style pppm 0.00001
kspace_modify slab 3.0
dielectric 1
#pair_modify shift yes #it is uncommented for shifted/truncated.

#Run/Print styles
fix 2 not_bottom langevin 1 1 6.99 904297 #not bottom is about 99% of total number of atoms

thermo_style custom step temp pe evdwl epair
thermo 10

run 100

Output:
With pair_modify yes
Step Temp PotEng E_vdwl E_pair
10285000 0.99044634 12.430849 0.0092685482 -0.23489335
10285010 0.99945363 12.41542 0.0070048401 -0.23700893
10285020 1.0001562 12.415137 0.0088646568 -0.23517801
10285030 0.98889901 12.432286 0.0086868613 -0.23548373
10285040 0.99976159 12.421335 0.008177021 -0.23611516
10285050 1.0043689 12.416264 0.0083009996 -0.23617625
10285060 0.99181344 12.437931 0.0099223853 -0.23468554
10285070 1.0077589 12.413876 0.0090538319 -0.23559688
10285080 1.0015782 12.424655 0.0075691142 -0.23696984
10285090 1.0042804 12.420705 0.0097466878 -0.23470351
10285100 1.0026287 12.426541 0.0094367655 -0.2346721
Loop time of 9.14612 on 1 procs for 100 steps with 6700 atoms

With no pair_modify yes
Step Temp PotEng E_vdwl E_pair
10285000 0.99044634 12.430849 0.0092685482 -0.23489335
10285010 0.99945363 12.41542 0.0070048401 -0.23700893
10285020 1.0001562 12.415137 0.0088646568 -0.23517801
10285030 0.98889901 12.432286 0.0086868613 -0.23548373
10285040 0.99976159 12.421335 0.008177021 -0.23611516
10285050 1.0043689 12.416264 0.0083009996 -0.23617625
10285060 0.99181344 12.437931 0.0099223853 -0.23468554
10285070 1.0077589 12.413876 0.0090538319 -0.23559688
10285080 1.0015782 12.424655 0.0075691142 -0.23696984
10285090 1.0042804 12.420705 0.0097466878 -0.23470351
10285100 1.0026287 12.426541 0.0094367655 -0.2346721
Loop time of 9.14048 on 1 procs for 100 steps with 6700 atoms

Total # of neighbors = 715493
Ave neighs/atom = 106.79
Ave special neighs/atom = 1.1606
Neighbor list builds = 5

Any help will be greatly appreciated.

Actually, I figured the error. Sorry to trouble. I read it from a restart file – so the pair_modify shift yes is automatically there. Sorry about that.

One other related question: If I use truncated vs shifted-truncated for the same initial conditions (apart from statistical uncertainties), will the final trajectory be identical (again neglecting numerical errors that may creep in)? If that is the case, won’t the difference in internal energy/atom from vdwl interactions between the two systems be equal to ave_neigh/atom*offset (approximately)?