# Potential Energy curve for CO2 and a surface

I am trying to calculate the potential energy curve for the interaction between CO2 and a silica surface. I use 3 systems and run single point calculations with LAMMPS and substract the values like E = (E_slab/co2 - E_loneco2 - E_lone_silica)

But the energy profile shows no going to zero at longer distances of CO2 away from the surface, can anyone Identify the reason?

my lammps input script for the surface+silica(total system) is as follows:

# LAMMPS Molecular Dynamics Simulation - single point calculation (crystoballite and Na+)

units           real
atom_style      full

boundary        p p p

pair_style     lj/cut/coul/long 7.00 7.00

pair_modify    mix arithmetic
pair_modify    tail yes
kspace_style   pppm 1.0e-4

bond_style     harmonic
angle_style    harmonic

#                          epsilon    sigma
pair_coeff      1       1  0.00010  2.70000  #SI-SI (bulk/silanol)
pair_coeff      2       2  0.18520  3.30000  #O-O (silanol)
pair_coeff      3       3  0.00000  0.00000  #H-H (silanol)
pair_coeff      4       4  0.46680  2.70000  #O-O (bulk)
pair_coeff      5       5  0.05589  2.75700  #C-C (CO2)
pair_coeff      6       6  0.15998  3.03300  #O-O (CO2)
pair_coeff      1       5  0.00236  2.72850  #Si_C(bulk_CO2)
pair_coeff      1       6  0.00399  2.86650
pair_coeff      2       5  0.21793  3.02850
pair_coeff      2       6  0.17212  3.16650
pair_coeff      3       5  0.00000  0.00000
pair_coeff      3       6  0.00000  0.00000
pair_coeff      4       5  0.16152  2.72850
pair_coeff      4       6  0.27327  2.86650
pair_coeff      5       6  0.09455  2.89500

bond_coeff      1 312.3295  1.66300    # Si-O    (surface) (SHAKE)
bond_coeff      2 553.0000  0.95000    # O-H               (SHAKE)
bond_coeff      3 575.0000  1.49900    # C-O     (CO2)     (SHAKE)

angle_coeff     1 24.49933  118.500    # SI-O-H    (surface)
angle_coeff     2 147.7000  180.000    # O-C-O     (CO2)

special_bonds   lj/coul 0.0 0.0 0.5

info            coeffs out log
neighbor        2.0 bin
neigh_modify    every 2 delay 0 check yes

thermo_style    custom step temp etotal ke pe density vol lx ly lz

run             0


The most likely culprit is using tail yes. That adds a correction energy which theoretically accounts for the truncation of the LJ parameter, but is only designed to work for systems of near-uniform density which is very different from what you are doing.