Energy conservation for Coulomb interaction models

Hi all,

My name is Neil and I use LAMMPS to perform simulations of ionic liquid electrospray. Currently I am trying to understand how Coulomb interactions are implemented in LAMMPS for the NVE simulations. I have browsed through some older posts in which Dr. Alex Kohlmeyer has mentioned that by definition, since the Coulomb model uses cut-off, these simulations cannot be considered as true energy conserving. This appears to be the case for a separate Coulomb model I have coded to LAMMPS for my own research work. However, when I use the coul/cut with a somewhat large cut-off radii, the NVE simulations run very smoothly, meaning that I do not see any energy gain. I have looked through the coul_pair_coul.cpp file but could not find any damping or any other smoothening being used. I was impressed by how well coul/cut was conserving energy. Does anyone have any pointers on why this is the case?
Thank you.

I am using latest version of lammps (13 March 2018) on STAMPEDE2 machine. I am simulating a 30,000 atom ethylammonium nitrate droplet, with an approximate diameter of 82 A. My input file is as follows:

All-atom simulation of ethyl-ammonium-nitrate

Neil Mehta, 2018

---------- Initialize Simulation ---------------------

clear
units metal
dimension 3
boundary f f f
atom_style full
#comm_modify cutoff 40.00
processors 4 4 4

---------- Other optional params ---------------------

bond_style harmonic
angle_style harmonic
dihedral_style opls

---------- Read the atom data -----------------------

read_data EAN_AA.inp

---------- Define Interatomic Potential --------------

pair_style hybrid/overlay lj/cut 12.0000095 coul/cut 40.00000
pair_coeff 1 1 lj/cut 0.0073717 3.250000
pair_coeff 1 2 lj/cut 0.0045931 3.372684
pair_coeff 1 3 lj/cut 0.0045931 3.372684
pair_coeff 1 4 lj/cut 0.0030967 2.850438
pair_coeff 1 5 lj/cut 0.0030967 2.850438
pair_coeff 1 6 lj/cut 0.0030967 2.850438
pair_coeff 1 7 lj/cut 0.0073717 3.199609
pair_coeff 1 8 lj/cut 0.0081932 3.048770
pair_coeff 2 2 lj/cut 0.0028619 3.500000
pair_coeff 2 3 lj/cut 0.0028619 3.500000
pair_coeff 2 4 lj/cut 0.0019295 2.958039
pair_coeff 2 5 lj/cut 0.0019295 2.958039
pair_coeff 2 6 lj/cut 0.0019295 2.958039
pair_coeff 2 7 lj/cut 0.0045931 3.320391
pair_coeff 2 8 lj/cut 0.0051050 3.163858
pair_coeff 3 3 lj/cut 0.0028619 3.500000
pair_coeff 3 4 lj/cut 0.0019295 2.958039
pair_coeff 3 5 lj/cut 0.0019295 2.958039
pair_coeff 3 6 lj/cut 0.0019295 2.958039
pair_coeff 3 7 lj/cut 0.0045931 3.320391
pair_coeff 3 8 lj/cut 0.0051050 3.163858
pair_coeff 4 4 lj/cut 0.0013009 2.500000
pair_coeff 4 5 lj/cut 0.0013009 2.500000
pair_coeff 4 6 lj/cut 0.0013009 2.500000
pair_coeff 4 7 lj/cut 0.0030967 2.806243
pair_coeff 4 8 lj/cut 0.0034418 2.673948
pair_coeff 5 5 lj/cut 0.0013009 2.500000
pair_coeff 5 6 lj/cut 0.0013009 2.500000
pair_coeff 5 7 lj/cut 0.0030967 2.806243
pair_coeff 5 8 lj/cut 0.0034418 2.673948
pair_coeff 6 6 lj/cut 0.0013009 2.500000
pair_coeff 6 7 lj/cut 0.0030967 2.806243
pair_coeff 6 8 lj/cut 0.0034418 2.673948
pair_coeff 7 7 lj/cut 0.0073717 3.150000
pair_coeff 7 8 lj/cut 0.0081932 3.001499
pair_coeff 8 8 lj/cut 0.0091063 2.860000
pair_coeff * * coul/cut
#kspace_style pppm 1.0e-05
#kspace_style barhut 20 2 40

neighbor 2.0 bin
neigh_modify delay 1 check yes page 800000
neigh_modify one 80000
balance 1.1 x 0.375 0.5 0.625 y 0.375 0.5 0.625 z 0.4843 0.5 0.5156

---------- Define Run-time settings ------------------

group ionic id 1:30000
velocity all create 5.0 4928459 rot yes dist gaussian
run_style verlet

---------- Writing the output ------------------------

dump 1 all xyz 100 traj.xyz