Grand canonical Monte Carlo simulations of a simple model semiflexible equilibrium polymer system, consisting of hard sphere monomers reversibly self-assembling into chains of arbitrary length, have been performed using a novel sampling method to add or remove multiple monomers during a single MC move. Systems with two different persistence lengths and a range of bond association constants have been studied. We find first-order lyotropic phase transitions between isotropic and nematic phases near the concentrations predicted by a statistical thermodynamic theory, but with significantly narrower coexistence regions. A possible contribution to the discrepancy between theory and simulation is that the length distribution of chains in the nematic phase is bi-exponential, differing from the simple exponential distribution found in the isotropic phase and predicted from a mean-field treatment of the nematic. The additional short length-scale characterizing the distribution appears to arise from the lower orientational order of short chains. The dependence of this length-scale on chemical potential, bond association constant, and total monomer concentration has been examined. (C) 2004 American Institute of Physics.