Protein crystallization conditions are usually identified by empirical screening methods because of the complexity of the process, such as the existence of nonequilibrium phases and the different crystal forms that may result from changes in solution conditions. Here the crystallization of a model protein is studied using computer simulation. The model consists of spheres that have both an isotropic interaction of short range and anisotropic interactions between patch-antipatch pairs. The free energy of a protein crystal is calculated using expanded ensemble simulations of the Einstein crystal, and NpT-Monte Carlo simulations with histogram reweighting are used to determine the fluid-solid coexistence. The histogram reweighting method is also used to trace out the complete coexistence curve, including multiple crystal phases, with varying reduced temperature, which corresponds to changing solution conditions. At a patch-antipatch interaction strength five times that of the isotropic interaction, the protein molecules form a stable simple cubic structure near room temperature, whereas an orientationally disordered face-centered-cubic structure is favored at higher temperatures. The anisotropic attractions also lead to a weak first-order transition between orientationally disordered and ordered face-centered-cubic structures at low temperature, although this transition is metastable. A complete phase diagram, including a fluid phase, three solid phases, and two triple points, is found for the six-patch protein model. A 12-patch protein model, consistent with the face-centered-cubic structure, leads to greater thermodynamic stability of the ordered phase. Metastable liquid-liquid phase equilibria for isotropic models with varying attraction tails are also predicted from Gibbs ensemble simulations. (C) 2004 American Institute of Physics.