Przedstawiamy metodę umożliwiającą dokładne przewidywanie harmonicznych energii swobodnych Helmholtza dla kryształów molekularnych w warunkach wysokoprzepustowych (high-throughput). Osiągnięto to poprzez opracowanie wydajnego obliczeniowo podejścia, które wykorzystuje regresję procesów gaussowskich (Gaussian Process Regression) opartą na lokalnych otoczeniach atomowych. Koszt trenowania modelu z wykorzystaniem potencjałów ab initio został znacznie zredukowany dzięki zastosowaniu potencjału empirycznego do początkowej optymalizacji parametrów modelu oraz doboru zbiorów treningowych i walidacyjnych. Dopiero potem model został transferowany i trenowany z użyciem danych opartych na teorii funkcjonału gęstości (DFT), z uwzględnieniem korekcji dyspersyjnych. Nasze rozwiązanie zostało przetestowane na zbiorze 444 struktur krystalicznych węglowodorów, obejmującym 38 polimorfów oraz 406 struktur zmierzonych w różnych warunkach lub pochodnych tych polimorfów. Osiągnięto wysoką dokładność predykcji i wydajność, uzyskując średnie bezwzględne odchylenie poniżej 0,04 kJ/mol na atom w temperaturze 300 K, przy trenowaniu modelu na zaledwie 60 strukturach krystalicznych. Ponadto, wykazaliśmy skuteczność i dokładność opracowanego podejścia, obliczając rozszerzalność cieplną krystalicznych węglowodorów aromatycznych w ramach przybliżenia quasi-harmonicznego oraz przewidując, jak zmiany objętości sieci wpływają na stabilność różnych polimorfów.