For instance, consider a system with $p_x$ and $p_z$ orbitals at the vertices of a square (on xy-plane). A square by itself would have $D_4$ symmetry. However, because of the $p_x$ orbital; the $90^\circ$ rotation ($C_4$) and $270^\circ$ rotation ($C_4^{-1}$) are no longer symmetry operations. Now the rest of the $D_4$ members do preserve the symmetry, but do not form a group since the subset is not closed.
How does one go about formulating a symmetry group for cases like these?
Hmm, so I'm not entirely sure what notation you're using ($D_4$ sounded like Schönflies to me), but I'm sure you can figure out the equivalent.
In general I agree with Andrew's point that molecular symmetry is determined by atomic positions and not orbitals. But at least to me as of now, this seems to be a thought exercise more than anything (after all, no real atom has only two p-orbitals), so I'll play along. In Schönflies notation your system would be labelled $D_\mathrm{2h}$. The three $C_2$ axes and three mirror planes are illustrated below. I left out the identity operation $E$ and the centre of inversion $i$.
