Potential Energy curve for CO2 and a surface

Question

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

read_data      data.lammps

#                          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

Answer 1

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.

Please note that without us knowing your data file or the potential profile you are observing, it is impossible to rule out any additional issues.