How to run crystal simulation on GAMESS?
As you said, there are many packages capable of handling calculations of crystals, slabs, etc. with periodic boundary conditions. VASP is only one of dozens. (That's probably a separate question in its own right.)

GAMESS does not include treatment of periodic boundary conditions, so you would need to run a calculation on a "cluster."

Since you mention wanting to treat a surface reaction, you'd use a builder (e.g., Avogadro) to create a large cluster (e.g., Build->Supercell). So rather than one primitive unit cell, you create maybe a 2x2x2 or 3x3x3 slab. If needed, you'll terminate dangling bonds (e.g., adding hydrogens to cap oxygens in $\ce{TiO2}$ for example.. notice the terminal "uncapped" atoms on the back and right of the unit cell).



Finally, after running calculations in GAMESS, you'd need to check for convergence. So you'd want to run calculations with increasing system size and see if the properties converge.

Now, this isn't as strange as it sounds. In some cases ($\ce{MgO}$ and $\ce{TiO2}$ are good examples) the properties require more advanced calculation methods like multi-configuration QM, and this may be easier to do in a quantum chemical package with a cluster, rather than a solid-state crystal code.