In Eq. 10, λ is a coupling parameter controlling the linear interpolation between the physical ensemble, characterised by potential energy function (U ), and the ensemble of the Einstein crystal, described by harmonic potential energy function ( k 2 ∥r -r 0 ∥ 2 ), where r 0 can be any static configuration of the system with centre of mass shifted to zero, and k is the scalar spring constant. We set r 0 to the initial configuration, with the centre of mass removed. The spring constant was fixed to 6 × 10 3 kJ(mol) -1 (nm) -2 so that the Gaussian distributions associated with the harmonic potential have negligible periodic images along the smallest box length (0.6358nm). The ground truth FE of a given supercell (f ECM crys ) is then estimated as: 55 11) where µ i = ( N j=1 m j ) -1 m i are normalised masses m i of the atoms in the system. In this work, n λ = 20 separate simulations were ran with λ parameters λ 1 , ..., λ 20 taking on values 0, ..., 1, respectively. The intermediate values, arranged in ascending order between zero and one, are reported in SI Figure . The same set of lambda values were used in all systems. Since the λ 1 = 1 simulation samples the translationally invariant ensemble of the unperturbed system of interest, all other simulations were initialised from r 0 with CMMotionRemover functionality of openMM being active. This removes average momentum of the centre of mass during the simulations. Additionally, every 50 timesteps during these simulations, the centre of mass was actively removed from the positions (r) to make sure that all systems remained on the same 3(N -1) dimensional subspace. In Eq. 11, ∆f i was calculated n λ -1 = 19 times using BAR functionality of pymbar. This was done several times during the simulations to show how the FE estimates converge as a function of more data becoming available for the BAR calculations (SI Figure ). The set of ground truth FEs reported in this paper utilised all of the available data. The standard error of f ECM crys was estimated by directly adding up the n λ -1 intermediate standard errors provided by pymbar. This approach provides a lower-bound estimate of the uncertainty on f ECM crys , because this error estimate assumes complete decorrelation between the samples used to compute free energy differences between neighboring replicas. While underestimating the uncertainty on the ground-truth method, this approach is fast (i.e. does not require bootstrapping) and provides a relatively more ambitious target for the PGM-based methods to compare with. Finally, f 0 appearing in Eq. 11 describes the normalisation constant of f ECM crys , respectively taking into account the FE of the Einstein crystal, the FE difference associated with removing the centre of mass from the Einstein crystal, and the FE difference associated with removing the centre of mass from the unperturbed physical system (where V was set to the volume of the box). Using the ground truth method described above (ECM), a total of 1.84 µs of data was sampled, with 0.32 µs of this data (17.4%) sampling the physical ensembles (i.e., λ = 1). This percentage is > 5% because each physical simulation was four times longer than any other simulation of the same ECM run to collect sufficient MD data for later use in the PGM models. Collecting this data in the context of ECM ensured that the average potential energy in each system was the same across all FE methods being compared in this work. Notably, one major methodological discrepancy between the ground truth method and the PGMs can be linked to how the centre of mass was removed. In the EC calculations, the centre of mass of the supercells was removed based on the position of all atoms (including the hydrogens). Meanwhile, in the PGMs, the centre of mass being removed was only based on the oxygen atoms, a more practical choice given the fact that the oxygens are used to anchor Cartesian blocks(Eq. 7). This methodological difference does introduce small systematic discrepancies between the absolute free energy of supercells, that cancel out when considering free energy differences between polymorphs, as shown in the Results section.