For each pressure, the system has been pre-equilibrated with 1.5 ns classical MD in the NVT ensemble at 300 K (while keeping fixed the BLG atoms), using the OPLS-AA force field. The latter has been extensively shown to accurately describe the structural organization of bulk methanol. Table . The first column displays the pressure conditions. Columns two, three and four report the optimized c cell parameter and the theoretical and experimental methanol densities In the next step, we performed Born-Oppenheimer MD based on the DFT-PBE approximation 60 together with Grimme's D3-BJ correction for dispersion interactions. This DFT level has been shown to reproduce within a 4.3% error the lattice constants, bulk moduli, elastic constants and interlayer binding energy of both graphite and diamond in the 0-25 GPa range. We employed Goedecker-Teter-Hutter pseudopotentials and combined plane-waves (850 Ry cut-off) and DZVP-MOLOPT SR (short range) basis sets, as implemented in CP2K. The electronic wavefunctions have been converged at each step imposing a threshold for the energy difference between two SCF cycles of 10 -6 Ha. All ab initio MD simulations are carried out in the spin-restricted Kohn-Sham scheme, with a timestep of 0.5 fs and the Nosé-Hoover thermostat 66 using a time constant of 100 fs, a chain length set to 3, the 3rd order Yoshida integrator and a multiple time step set to 2. Each BLG/methanol system was simulated at 300 K for 20 ps: the first 5 ps are considered as equilibration, while structural analysis (density profile and radial distribution functions) was performed on the last 15 ps of each trajectory. Fig. in Supporting Information reports the distribution of stacking configurations explored during MD between 0.0001-10 GPa.