Identifying the solid phase is the easiest. you will see a crystal, an organized structure with six neighbors for each particle. (the Mermin–Wagner theorem does no apply here, because the interaction has a long enough range. see this article for simulations and snapshots of the solid phase). here is what the solid phase looks like:
there will be a hexatic phase (liquid crystal, I can't find snapshots).
The hard part is distinguishing the liquid and the gas. you should notice that for a high enough temperature there is no distinction between gas and liquid (above the critical point) so if you just look at the distance between particles, you can't say what's liquid and what is gas.
If you look at a temperature lower than the critical temperature, than depending on the pressure you apply, you can see gas, liquid, or drops of liquid within a gas. when I say drops of liquid within a gas, I mean that you will have clusters of particles and some free particles moving between them - see the picture in the liquid gas region.
if you are above the critical temperature there is no distinction between gas and liquid, we call this fluid (see the one-phase region in the picture).
If you want to distinguish between gas \ liquid \ fluid - there is a way:
what you can do is look at the density distribution. divide a snapshot of your simulation into boxes (you will need to play with the sizes, depending on the number of particles and density), like this picture:
ask what is the density in each box, and put the results in a histogram. do this for all the steps in your simulation, and average the results. see this article for more details.
if you have liquid mixed with gas, you will see two densities:
as you increase the temperature, the two densities will get closer, until they become one. this means you crossed the critical temperature, and you are now above the critical temperature.
one more note: finding the freezing point and the vaporization point can require different simulation techniques. as you get closer to the critical point, you will need much more steps to get the simulation equilibrated, so you will probably want to use the Gibbs ansamble to avoid this.