Born-Oppenheimer molecular dynamics (BOMD) simulation controlled by using the Nosé-Hoover thermostat was performed at a temperature of 500 K for 5 ps to identify the thermal stability of the designed structures. The Heyd-Scuseria-Ernzerhof screened hybrid functional (HSE06) was employed to calibrate the band structures as the GGA-PBE functional usually underestimates the band gaps. The accuracy of HSE06 in reproducing the band structure of 2D HT-polymers has been verified in previous studies. In order to assess the ability of photocatalytic water splitting of 2D binary polymers, the band edges were corrected by the vacuum level and compared with the water redox potentials. Carrier mobilities were calculated using deformation potential theory (DPT). Detailed computational settings are given in supporting information. 2D binary HT-polymers show a reduced D3h symmetry, analogous to that of h-BN. The optimized lattice parameters are found to be 17.47, 17.15 and 17.53 Å for CTPAB, OTPAB and MTPAB, respectively, which are very similar and inbetween those of the unary B and N centered HT-polymers. More structural details including pore size and distance of two in-plane neighboring center atoms are also exhibited in Table . BOMD simulations (Figure ) show that 2D binary HT-polymers remain intact at 500 K after 5 ps and recover the initial structure after subsequent geometry optimization, indicating the high thermal stability of the designed structures. The integrity of the alternating arrangement of B and N heteroatoms is confirmed by calculating the relative formation energy with respect to their monomers: