. The resulting equations were numerically propagated using the split-operator method 37 on the grid R = (3 a.u. -100 a.u.) with the total of 1024 points. In all simulations we explicitly check whether reaches the end of the grid. Such an occurrence would correspond to a spurious appearance of the part of the wavepacket at 3 a.u. (since the split-operator method employs Fast Fourier Transform with periodic boundary conditions) and in turn unphysically distort the nuclear dynamics. We performed several numerical experiments and found that for our choice of parameters, if the R-grid is terminated at 100 a.u., simulations can safely be performed up to 500 fs. The dipole moment, , of each molecule is evaluated at every time step and the macroscopic polarization is computed following the mean field approximation 38 (13) where is the number density of molecules in the slab. The polarization is coupled back to Maxwell's equations ( ) via polarization current (last term in Eq. (9b). The resulting system comprised of Eqs. ( ), (10), and ( ) is solved numerically as described above.