# Electric Field and Temperature Effects on the Ab Initio Spectroscopy of Liquid Methanol
## 2 Materials and Methods  
We used the software package CP2K [51; 52], based on the Born-Oppenheimer approach, to perform ab initio molecular dynamics (AIMD) simulations of samples of neat liquid methanol under the action of static and homogeneous electric fields applied along a given direction (corresponding to the \(z\)-axis). Moreover, in the zero-field case, several AIMD simulations have also been carried out at different temperatures. The implementation of an external field in simulation algorithms based on density functional theory (DFT) can be achieved by exploiting the Modern Theory of Polarization and [PERSON]'s phases [53; 54; 55], whereas when periodic boundary conditions apply--as in the current investigation--implementation of static fields is executed owing to the seminal work by [PERSON] and [PERSON] [56]. The interested reader may refer to [34] for a vast presentation of the most used implementations of electric fields in atomistic simulation codes. Neat methanol samples contained 50 CH\({}_{3}\)OH molecules (i.e., 300 atoms) arranged in a cubic cell with edges equal to 14.99 A, so as to reproduce the liquid methanol experimental density of 0.79 g\(\cdot\)cm\({}^{-3}\) at room temperature. As anticipated above, in order to minimize nonphysical surface effects, such a structure was replicated in all spatial directions by adopting periodic boundary conditions. The intensity of the electric field was gradually increased with a step increment of 0.05 V/A from zero up to a maximum of 0.30 V/A, the latter value being close to the electric field threshold for molecular dissociation at 300 K [49; 50]. Beyond this field intensity, indeed, a plethora of electric-field-induced chemical reactions takes place in liquid methanol [57; 58]. Moreover, when evaluating pure temperature-induced effects, four different AIMD simulations were executed at the nominal temperatures of 240, 270, 300, and 350 K. In these zero-field cases we performed dynamics of 50 ps, whereas for each other value of the field intensity, we ran dynamics of 20 ps, thus accumulating a global simulation time equal to 320 ps, whereas a time-step of 0.5 fs has been chosen. Wavefunctions of the atomic species were expanded in triple-zeta valence plus polarization (TZVP) basis sets with Goedecker-Teter-Hutter pseudopotentials using the GPW method [59]. A plane-wave cutoff of 400 Ry was imposed. Exchange and correlation effects were treated with the Becke-Lee-Yang-Parr (BLYP) [60; 61] density functional. Moreover, in order to take into account dispersion interactions, we employed the dispersion-corrected version of BLYP (i.e., BLYP-D3) [62; 63]. Adoption of the BLYP-D3 functional has been indicated by the widespread evidence that such a functional, when dispersion corrections are taken into account, offers some of the best adherence to the experimental results related to H-bonded systems like water among the standard generalized gradient approximation (GGA) functionals [64; 65]. It is well-known indeed that neglecting dispersion corrections leads to a severely over-structured liquid (see, e.g., [66] and references therein). Although some recent findings have magnified the importance of including nuclear quantum effects (NQEs) in the presence of electric fields [33], the required computational cost for the wide palette of simulations here presented renders impracticable the adoption of sophisticated thermostats mimicking NQEs. More importantly, vibrational spectra determined from path-integral AIMD exhibit many critical issues due to the presence of fictitious vibrations induced by thermostats for reproducing NQEs [67; 68]. This way, the dynamics of nuclei were simulated classically within a constant number, volume, and temperature (NVT) ensemble using the Verlet algorithm, whereas the canonical sampling was executed by employing a canonical-sampling-through-velocity-rescaling thermostat [69] set with a time constant equal to 10 fs.