We don’t use VASP’s elastic constant routine. By using VASP’s built-in elastic routine, you’re doing something different than what we perform and asking why you’re getting a different answer. As Joey mentioned, our workflow is in atomate, and you’re welcome to run it directly and if you’re still getting different results we could help. You can also try manually running our procedure as follows:

1.) Relax the structure to an energy convergence of ~ 1E-6 using a K-point density of 7,000 per reciprocal Angstrom^3 and turn off symmetry so VASP can’t mess with the K-mesh.

This has to be done in a separate calculation. Are you doing this?

2.) Generate 4 strained structures corresponding to -1%, -0.5%, 0.5%, 1% for each of the 6 unique deformation directions. I believe VASP only uses 1 per direction.

3.) Run static calculations for each of these new deformed structures using an ENCUT of 700, k-point density of 7,000 per reciprocal angstrom^3, turning off symmetry, using ADDGRID=True to get an even more precise stress tensor. In addition, the K-point mesh for these calculations is for the symmetry of the deformed structure. VASP uses the k-point mesh of the original structure, as far as I can tell.

4.) We then fit these back to the elastic tensor using an over-determined matrix inverse, which lets us separate out any potential non-linearities.

Our procedure is excessive when a material is hard and insulating, such as mp-48 where you are likely to get a good elastic tensor with a simple procedure as what is built into VASP. As it gets softer and more semiconducting, the effect of all these differences gets magnified. Especially for something that is potentially soft. For the case of PbGeO3 in particular, I would be worried about the strains that VASP uses for generating the deformations, whether you have a high enough ENCUT + ADDGRID to get accurate stress tensors and whether you’re relaxing sufficiently to compute an elastic tensor in the first place over the other INCAR parameters.