I have checked my calculations, the k points and cutoff energies have been carefully tested. I wonder if the error comes from the basis sets we used (DZP)? Is there anyone who used to perform similar calculations? Is that a common error in SIESTA? I urgently need your help, thank you!
Try to reproduce other properties first (lattice constant, bulk modulus, etc.) and see if they are reproduced well - that would tell you about the quality, then you will know whether you need to refine the basis or maybe even pseudopotential.
It is also important how much values and what strain do you impose to use for fitting. Strain should not exceed 1%.
Bulk modulus is not well reproduced in DFT, for elastic constants the error is even higher. There is nothing so special in SIESTA to prevent you from obtaining values of accuracy comparable to other packages with pseudopotential approach.