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