To estimate electromechanical response of the FE doped h-BN, we calculate piezoelectric stress constant e 33 from the change of the polarization along the cdirection, normal to the h-BN plane, as ∆P = P sp 3 (η 3 ) -P sp 3 (η 3 = 0), in response to a uniaxial strain along the same direction (η 3 ). Figure Table . Ionic (e ion 33 ) and electronic (e elc 33 ) parts of total e33obtained from the self-consistent response to a finite electric field ε (LCALCEPS-tag in VASP) are in C/m 2 unit. The piezoelectric constants obtained from DFPT calculations using the optimized structures considering van der Waals interaction are shown in the parenthesis. Elastic constant C33 and longitudinal piezoelectric strain constant d33 are in GPa and pC/N unit, respectively. Note that 1 C/N equals to 1 m/V. to η 3 . Using this slope method, our calculated e 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.282 C/m 2 , -0.912 C/m 2 , -0.306 C/m 2 , and -1.169 C/m 2 , respectively. Interestingly, all the FE structures show negative sign for e 33 as P sp 3 decreases in response to elongation of c lattice parameter (i.e. positive η 3 ).