Interfacial operator approach to computing modes of surface plasmon polaritons for periodic structures