To understand the origin of the negative sign, we decompose e 33 into two terms: e ion 33 where the piezoelectric response is only due to ionic displacement in response to η 3 , and e elc 33 which comes only from redistribution of electron density in response to η 3 . We calculate e elc 33 from the slope of the linear relation between ∆P and η 3 (see Figure ), where atoms are kept fixed at their zero strain optimized positions and only the electron density is optimized in response to η 3 . Our calculated e elc 33 (from which e ion 33 = e 33 -e elc 33 ) for Al 0.0625 B 0.9375 N, Al 0.1250 B 0.8750 N, Ga 0.0625 B 0.9375 N, and Ga 0.1250 B 0.8750 N are -0.078 (-0.154) C/m 2 , -0.251 (-0.661) C/m 2 , -0.169 (-0.137) C/m 2 , and -0.395 (-0.774) C/m 2 , respectively. We notice that the e ion 33 are much larger than the e elc 33 . However, both e elc 33 and e ion 33 have negative sign, whcih is different from other piezoelectrics that show negative stress constants. For example, FE HfO 2 (orhtorhombic phase with Pca2 1 space group) shows negative e 33 (-1.53 C/m 2 ) due to large negative e ion 33 (-2.16 C/m 2 ) but positive small e elc 33 (0.63 C/m 2 ) . On the other hand, van der Waals layered FE CuInP 2 S 6 exhibits a small negative e 33 ( -9.7 µC/cm 2 ) because of a large negative e elc 33 (-10.1 µC/cm 2 ) and a small positive e ion 33 (0.4 µC/cm 2 ) . Overall, our calculated e 33 are an order of magnitude larger than that of FE CuInP 2 S 6 .