There is a slightly confusing point about that tutorial - the part where you fit and later fix the c-value, uses FHI pseudopotentials/basis sets, while the first part (where you let c adjust itself self-consistently) is the part using the HGH basis sets. If you run it with just
exchange_correlation = MGGA.TB09LDA
and the HGH basis sets, the results look fine.