It all depends what information you're trying to get from the simulation. Do you want the iso-frequency contours to study bandgaps or possible regions of self-collimation? Then use something that implements the plane-wave expansion method (or implement it yourself in your language of choice should you ne inclined). Do you want to know what happens with an incident wave in a steady-state regime? Use a package that implements the finite-difference frequency-domain method. It is also important to note that almost any method/package can do whatever you need it to, they might just not be doing it efficiently. As other people have mentioned, for FDTD, MEEP and Lumerical work really well in 3D cases even.