temperature dependent bulk modulus from QHA-fit

[freshwind]

I have a problem with deriving the temperature dependent bulk modulus within the quasi-harmonic approximation. I did calculations on LiCoO2 with strict convergence cirteria for relaxations and force calculations. However, according to the quasi-harmonic fit the bulk modulus increases with temperature due to the increasing curvature (second derivative of the energy with respect to the volume) of the free energy with increasing temperature. I already tried increasing supercel size and the kpoints per reciprocal atom, to no effect.This is somewhat contradictory to physical reasoning and I don't know of materials with that behaviour. I surmise this is a problem of the PAW potentials since they effectively determine the volume dependence of the total energy and the vibrational free energy.
Can someone offer a constructive solution?