Interaction between similarly charged surfaces can be attractive at high electrostatic coupling constants = lBZ2/μGC, where lB is the Bjerrum length, μGC the Gouy-Chapman length, and Z the valency of counterions. While this effect has been studied previously in detail, as a function of surface charge density and valency of the pointlike counterions, much less is known about the effect of counterion size. We apply the Wang-Landau sampling Monte Carlo (MC) simulation method to compute the free energy F as a function of the scaled distance between the plates D̃ =D/ μGC for a range of and scaled counterion radii R̃ =R/ μGC. We find that for large and small ion radius, there is a global equilibrium distance D̃ = D̃ eq =2 (1+ R̃), correctly giving the expected value at the point counterion limit. With increasing R̃ the global minimum in F(D̃) changes to a metastable state and finally this minimum vanishes when R̃ reaches a critical value, which depends on . We present a state diagram indicating approximate boundaries between these three regimes. The Wang-Landau MC method, as it is applied here, offers a possibility to study a wide spectrum of extended problems, which cannot be treated by the use of contact value theorem. © 2010 American Institute of Physics.