I'll try to outline an approach in plane-wave DFT. The main idea is to build up the system step-wise and reuse previous results as the number of atoms and therefore the computational burden will grow quickly.
Lets start with the simulation using a bulk material with $N$ atoms and assume the surface direction is orientated facing the z-axis. Make sure this cell is well relaxed.
Build a super cell (SC1) by repeating the cell along the z-direction $(2n+1)$-times that will have $(2n+1)N$ atoms.
An initial guess for the charge of the SC1 can be done by periodically repeating the charge from the first simulation. Eventually you need to interpolate this density such that it fits to the grid of the DFT simulation for SC1. In this simulations that atoms do not need to be relaxed as it's basically still bulk.
Now include a vacuum in the z-direction. Again an initial guess can be made from the charge density of the last simulation but the grid must be carefully adapted. Here, I would relax the first top and bottom layer and freeze all other layers, but only in z-direction. If they move a lot, make sure the vacuum included is sufficient. If so restart with relaxing 2 upper and bottom layer.
Once this is done make sure that the charge density in the middle layer is still bulk-like. If not, restart from step 2 with more layers.
As we want to study the interaction with a molecule we need a wide enough cell in x and y-direction to avoid the molecule interact with its periodic replica. Thus we build a second super cell (SC2) by repeating the structure of the last simulation $n_x$-times in $x-$ and and $n_y$ times in $y$-direction. Again, reuse the charge density from the last simulation as we now have $(2n+1)n_x n_y N$ atoms.
Finally, place your molecule close enough to one of the surface and relax only the atoms close to it.
Some practical hints in VASP,
The python function resample from scipy.signal is quite useful to interpolate the charge density.
Be aware of augmentation charges in VASP when extending your charge.
NGX, NGY and NGZ are the parameters that specify the real space grids.
To analyse what is happening I would consider the potential that can be obtained with the Tag LVTOT = T and the differences with and without the molecule.
The differences in charge density will tell further details about the interaction. So I would eventually simulate the molecule alone in SC2 and the compare to the situation including the surface.