surface energy

Hi Dear lammps users
I am calculating size dependent surface energy of copper.
i use this formula:
surface energy=(E ,slab)_(N*E,bulk)/2A
where E bulk is the bulk cohesive energy that i calculated with this code:

Find minimum energy fcc configuration

Mark Tschopp, 2010

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

clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array

---------- Create Atoms ---------------------

lattice fcc 3.61
region box block -5 5 -5 5 -5 5 units lattice
create_box 1 box

lattice fcc 3.61 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 box
mass 1 63.546
replicate 1 1 1

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

pair_style eam
pair_coeff * * cu_u3.eam
neighbor 2.0 bin
neigh_modify delay 10 check yes

---------- Define Settings ---------------------

compute eng all pe/atom
compute eatoms all reduce sum c_eng

---------- Run Minimization ---------------------

reset_timestep 0
fix 1 all box/relax iso 0.0 vmax 0.001
thermo 10
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms
min_style cg
minimize 1e-25 1e-25 5000 10000

variable natoms equal “count(all)”
variable teng equal “c_eatoms”
variable length equal “lx”
variable ecoh equal “v_teng/v_natoms”

print “Total energy (eV) = {teng};" print "Number of atoms = {natoms};”
print “Lattice constant (Angstoms) = {length};" print "Cohesive energy (eV) = {ecoh};”

print “All done!”

now i want to calculate E slab.
so i should use 2 dimensional boundary conditions.
and i don’t know how to calculate the energy of slab.
I will be grateful if somebody help me
Thanks