We have studied the thermodynamic properties of hen egg white lysozyme crystals using a novel simulation method combining atomistic Monte Carlo simulation to calculate van der Waals interactions and the boundary element method to solve the Poisson-Boltzmann equation for the electrostatic interactions. For computational simplicity, we treat the protein as a rigid body, using the crystallographic coordinates of all non-hydrogen atoms of the protein to describe the detailed shape. NVT Monte Carlo simulations are carried out for tetragonal and orthorhombic crystals to obtain the van der Waals energy, incorporating an implicit solvation effect. For crystal phases, an optimally linearized Poisson-Boltzmann equation is used to include the effect of the Donnan equilibrium of the salt ions. The Helmholtz energy is obtained from expanded ensemble Monte Carlo simulations. By using the force field parameters that had previously been tuned for the solution properties, reasonable agreement with experiment is found for the crystallization energy of the tetragonal form. The prediction of the entropy is also reasonable with a slight underestimation suggesting the release of a few water molecules per protein during the crystallization. However, the predictions of the properties of the orthorhombic crystal are poor, probably due to differences in the solvation structure as indicated by experiments, and also as a result of the approximate force field used.