Photosystem II (PSII) and its associated light-harvesting complex II (LHCII) are highly concentrated in the stacked grana regions of photosynthetic thylakoid membranes. PSII-LHCII supercomplexes can be arranged in disordered packings, ordered arrays, or mixtures thereof. The physical driving forces underlying array formation are unknown, complicating attempts to determine a possible functional role for arrays in regulating light harvesting or energy conversion efficiency. Here, we introduce a coarse-grained model of protein interactions in coupled photosynthetic membranes, focusing on just two particle types that feature simple shapes and potential energies motivated by structural studies. Reporting on computer simulations of the model’s equilibrium fluctuations, we demonstrate its success in reproducing diverse structural features observed in experiments, including extended PSII-LHCII arrays. Free energy calculations reveal that the appearance of arrays marks a phase transition from the disordered fluid state to a system-spanning crystal. The predicted region of fluid-crystal coexistence is broad, encompassing much of the physiologically relevant parameter regime; we propose experiments that could test this prediction. Our results suggest that grana membranes lie at or near phase coexistence, conferring significant structural and functional flexibility to this densely packed membrane protein system.