For what concerns the parallel replica CG simulations, the equilibration was run for 2000 steps (200 ps) while the differentiable MD for 200 steps (20 ps). Since we use 40 replicas for each optimization, the gradients are effectively computed from 40 × 20 ps = 800 ps long simulations. In both equilibration and production, we used a time step of 0.02 ps for the interatomic forces calculation and a time step of 0.1 ps for the intermolecular and electrostatic forces. The systems were kept at constant temperature with the velocity rescale thermostat and at constant pressure with the semi-isotropic Berendsen barostat . The coupling constant was set to 0.1 ps for both. We used a 20 × 20 × 30 grid for the particle mesh calculations, compressibility κ = 0.05 kJ -1 mol, particle spread σ = 0.5 nm, ρ 0 = 8.33 nm -3 and a = 9.21 nm -3 . The bonded parameters used for the hPF model are reported in Table of the Supporting Information.