All of the DFT simulations are performed using the open source code SIESTA . The core electrons are replaced by norm conserving pseudopotentials following a Troullier and Martins scheme ; the valence electrons are represented by an extended numerical atomic basis set of the polarized double-zeta (DZP) type. The pseudopotential used each P atom was 3s 2 3p 3 3d 0 4f 0 accordingly. For the exchange-correlation part of the electron energy, we used the Generalized Gradient Approximation (GGA) with the PBE functional . The energy shift parameter used for the basis was 50 meV and the mesh cutoff was 500.0 Ry. For convergence of the self-consistent fields, we used an energy threshold of 1.0 × 10 -5 eV per supercell and a force-tolerance of 0.03 eV/ Å. The relaxation of atoms and supercell is carried out by using the conjugate gradient optimizer with a maximum displacement of 0.05 Å at each relaxation step and the maximum stress tolerance for each component of the stress tensor is set to be 0.1 GPa. The deformation simulations are performed by applying uniaxial strain at a strain rate of 0.005 along the loading direction up to ϵ = 0.5; at each of the quasi-static loading steps, the atoms and the supercell along the lateral directions are allowed to relax, which may lead to minor discrepencies between calculated macroscopic strain and intended strain. Throughout our results and discussions, we will refer to ϵ according to strain found by directly calculating the change in layer size or bond lengths.