Recent experiments have shown that exfoliated few-layer CrI3, a prototypical van der Waals magnet, undergoes a phase transition from the high-temperature monoclinic structure to the low-temperature rhombohedral structure under pressure. To understand how the magnetism of CrI3 responds to these structural changes, we perform ab initio density functional theory simulations on bilayer CrI3. We simulate the interlayer lateral shift-dependent potential energy surface of bilayer CrI3 to examine the stability and magnetism as a function of external pressure. Using the hybrid PBE0 functional, we are able to give qualitatively correct exchange coupling energies, without using an on-site Coulomb interaction correction. Thus, we avoid using tunable parameters. The results show that under external pressure, the monoclinic crystal structure is destabilized in comparison with the rhombohedral structure, in agreement with the observed phase transition in few-layer CrI3 devices under pressure. We also look into the microscopic origins of the interlayer exchange coupling. We identify the competing orbital pathways that favor ferromagnetic and antiferromagnetic kinetic exchange, respectively, which are consistent with previous reports. This study opens a new direction of using hybrid functionals with Gaussian orbitals and a cluster-based approach for obtaining Heisenberg J values to accurately simulate the magnetic properties of 2D materials.