**4. Computational methodology**

First-principles calculation was performed by using VASP code [34]. Interactions between ion and electron were treated by projector augmented wave (PAW) method [45]. PBEsol functional [46] was used to approximate exchanges and correlate interactions of electrons, which can be used to reproduce the lattice constants of various materials [45]. Precise calculation on the lattice constant is essential to predict piezoelectric properties because they depend on volume of unit cell Ω as shown in Eq. (1). The kinetic energy cutoff for plane waves was set at 500 eV, and the k-point mesh was set at ~0.03/Å intervals to obtain the converged total energy at less than 0.1 meV/atom. Before calculating the piezoelectric constants, atomic positions and cell parameters were optimized until the forces on each atom and cell converged at below 5 × 10−4 eV/Å.

Since VASP does not directly calculate Eq. (12), we added routine to calculate displacementresponse internal-strain tensor *Γnj* and decompose piezoelectric constants into each atom. The sum of the decomposed piezoelectric constants was confirmed to accurately reproduce the total piezoelectric constants. Careful convergence tests with a higher energy cutoff and denser *k*-point mesh showed that the numerical accuracy of the calculated *Γnj* was less than 0.01. It was confirmed that this error does not influence our discussion and conclusion. Moreover, it was confirmed that the values of *Γnj* obtained by the DFPT method were consistent with those calculated by the direct method, in which the strain-displacement relation of each ion was explicitly calculated.

On the basis of cubic *Pm3m* phase, lattice instability analysis was performed by phonon calculation utilizing phonopy code [47]. Force constant matrix shown in Eq. (2) was constructed by DFPT calculation implemented in VASP code combined with supercell approach. Supercell was constructed by using unit cell so that orthogonal three axes of the supercell exceed 10 Å. Note that although supercell is not required in DFPT approach, the present VASP code implements perturbation at the zone center.
