I have performed the calculations for Tetragonal BaTiO3 using QE.

You can find the results in the paper. (J. Phys.: Condens. Matter 32 (2020) 475601 (8pp))

"https://iopscience.iop.org/article/10.1088/1361-648X/abaa81/meta"

In case of Tetragonal BaTiO3 I choose z-direction as gdir , using the flags in nscf calculation

lberry = .true.

gdir = 3

nppstr = 24

also increase the K-points along the gdir direction.

For scf calculations I used

K_POINTS automatic

8 8 8 0 0 0

for nscf calculations I used

K_POINTS automatic

8 8 24 0 0 0

nppstr = number of k-points to be calculated along each symmetry-reduced string.