While doing this, I used a k-mesh of 8x8x8 and ecut = 120 Ry. I used the realxed lattice parameter 7.2794 bohr(7.288 bohr in the paper). With this, I got the lowest optical mode frequency as 42.46 cm^(-1).
The flexoelectric coefficients I got was (in nC/m),
xx,xx = -212.529 (-182.841), xx,yy = -184.381 (-158.785), xy,xy = -22.502 (-19.392). The values in the bracket are the coefficients taken from the above-paper. Is this difference in values acceptable? Does this have anything to do with the imposed boundary conditions?
I am attaching the input files of the calculations here. flexo.abi (2.6 KB) flexo.abo (417.7 KB) tlw_3.abo (23.1 KB)
I have compared your output files with the ones we produced four our paper, and the only relevant difference appears to be in the lattice-mediated contribution to the flexoelectric tensor. In particular, when the clamped-ion flexoelectric force-response tensor (which looks pretty similar in both datasets) is converted to a flexoelectric internal strain tensor by applying the pseudoinverse of the interatomic force constants. I think that such a differences arises because your calculation, due to using different pseudos, yields a lower energy for the soft polar phonon than ours. This is also the reason why your flexoelectric coefficients are larger in magnitude.
Apart from that, I recall that linear-response calculations in insulators should be done while using only valence bands. Your choice nband=23, with 3 empty bands, might be harmless in cubic STO but problematic in other systems.
Thanks for the clarification. In non-cubic systems, should I then go with
nbands = number of occupied bands +1 for the calculations, without specifying nbdbuf variable?