Here, we report the construction of a model Hamiltonian based on molecular pair interactions fitted to first principles calculated energies. Applying the model Hamiltonian, we perform Monte Carlo simulations to gain insight into the orderdisorder phase transition as well as nanoscale effects and defective phases that are likely to be involved in the piezoelectric performance. Particularly, we show the importance of including contributions from vibrational entropy to model phase transition temperatures in quantitative agreement with experiment. Finally, we show that the easy development of disorder in the ab-plane results in formation of isolated and extended defects, which are likely to be responsible for the exceptional piezoelectric response.