Estrogen receptor alpha (ERalpha) and beta (ERbeta) are ligand activated transcription factors that have different physiological functions and differential expression in certain tissues. The ligand binding domain of ERbeta shares 58% sequence identity with that of ERalpha. However, in the binding pocket there are only two relative residue substitutions. This high similarity at the active site is a great challenge for designing selective estrogen receptor modulators. ERbeta is shown to be related to several diseases. To understand the molecular basis of ERbeta selectivity, molecular dynamics simulations were carried out for both ERalpha and ERbeta complexes. Our simulations revealed the conformational changes at the active site of the ERs and the difference of affinity with ligand. The electrostatic repulsion between the Sdelta atom of ERalpha Met421 and the acetonitrile group nitrogen atom of the ligand led to unfavorable binding. The repulsion resulted in the conformational change of the side chain of ERalpha Met421, which changed the conformation of both Leu346 and Phe425. These residues changes expanded the volume of binding cavity, which led to unstable binding of the ligand. In addition, ERbeta Met336 formed more hydrophobic contacts with the ligand relative to corresponding residue ERalpha Leu384. Furthermore, the binding free energy analysis was shown to be correlated with the previous results determined by experiment. At last, free energy decomposition evidently indicated the contributions of key residues. The present results could help explain the mechanism of ERbeta selectivity and may be considered in the design of subtype-selective ligands.