All calculations utilised the revised version of the Peintinger-Oliverira-Bredow split-valence triple-ζ basis set (pob-TZVP-rev2) and the Becke-3-Lee-Yang-Parr (B3LYP) hybrid density functional. The B3LYP density functional was supplemented with Grimme's noncovalent dispersion correction (D3) and the Becke-Johnson damping correction with three-body Axilrod-Teller-Muto repulsion contributions (program keyword "ABC"). In all MgSO 4 heptahydrate calculations, 125 k-points were used in the irreducible Brillouin zone (keyword SHRINK=9) and 99 radial point and 1454 angular points were used for the pruned integration grid. The overlap-based truncation criteria for the bielectronic integrals (Coulomb and exchange) (program keyword "TOLINTEG") were set to 10 -12 , 10 -12 , 10 -12 , 10 -20 and 10 -40 for all calculations and the maximum order of multipolar expansion was set to 6 (program keyword "POLEORDR"). The starting structure for MgSO 4 heptahydrate was published by Ferraris, Jones, and Yerkess in 1973 and the initial ionic charges were explicitly set to Mg 2+ and (SO 4 ) 2-. In the geometry optimisation, the lattice dimensions and atomic positions were allowed to fully optimise within the P 2 1 2 1 2 1 space group, and the energy convergence was set to ∆E < 10 -8 E h . The optimised structure was used to calculate the vibrational frequency analysis and the energy convergence was set to ∆E < 10 -10 E h .