## 2 Computational Methods  
MD-MTD simulations allow for efficient sampling of free energy surfaces of two interacting molecules without the formation and breakage of a covalent bond in an accelerated time scale of less than one microsecond (\(\upmu\)s) by filling the potential wells with repulsive bias potentials [37; 38; 39; 40]. As shown in Figure 2, the bias potential is added to the local minima to facilitate barrier crossing. Once one local minimum is filled close enough to the barrier, the system crosses the barrier and moves to the neighboring minimum. When all the local minima are filled, the FES becomes flat and the system can move back and forth without any barrier. The FES is subsequently reconstructed based on the amount of bias potentials added. This method assumes that several collective variables (CVs) can be used to describe the process, distinguishing different states of the interacting system. Many variables including the number of hydrogen bonds formed, the coordination number, bond distance, bond angle, and dihedral angles, can be used as CVs. To describe the interacting system consisting of one polymer chain from the catalyst and one cellulose substrate accurately, the number of formed hydrogen bonds between the polymer chain and cellulose substrate was used as the CV, as shown in Equation (1). Because the catalytic polymer chain can be both the donor and the acceptor of a hydrogen bond, two CVs were used. One CV was the hydrogen bonds formed between one polymer chain (donor, group i), and one cellulose chain (acceptor, group j), and the other CV had the opposite setting of the donor and the acceptor. \[CV_{hydrogen\ bonds}\ =\ \sum_{ij}\frac{1-\left(\frac{d_{ij}}{r_{0}}\right)^{n}}{1 -\left(\frac{d_{ij}}{r_{0}}\right)^{m}} \tag{1}\]  
where \(d_{ij}\) is the distance between atom \(i\) and atom \(j\); \(r_{0}\) is the hydrogen bond cutoff distance. \(m\) and \(n\) are high-power integers used to determine the formation of hydrogen bond. The values \(m=12\), \(n=6\), and \(r_{0}=2.5\) A are typically chosen for calculating the hydrogen bonds. It is found that polymer chain length plays a critical role in the catalytic activity based on our previous experiment results [34]. MD simulations were conducted for two different chain lengths with degree of polymerization (DP) of 10 and 20 for PSSA, PIL, and cellulose chains. Earlier studies [41; 42; 43; 44] showed that the properties of ionic liquids (ILs) using an OPLS (optimized potentials for liquid simulations) force field from MD simulations agreed well with experimental results. The OPLS force field [45], if available, was used for PSSA and PIL chains for the catalysts. Additional force fields parameters of PSSA, PIL, and [EMIJM]Cl were parameterized based on quantum mechanics (QM) calculations using SCAN in Gaussian09 [46]. The monomer units in PSSA and PIL chains were terminated by two methyl groups to mimic the polymeric environment. The scan ranges of bond length, bond angle and dihedrals were 0.2 A, 10\({}^{\circ}\), and 360\({}^{\circ}\) respectively. Geometries for the neutral species were optimized at the MP2/6-31G**+ level in gas phase. The single point energy was subsequently calculated with the same method and the same basis set for each scanned structure. For PIL and ionized PSSA, the force field parameters were determined at the same level, and the same basis set  
Figure 2: Filling a one-dimensional model potential with bias potentials, starting in the right local minimum. After filling up the local minimum, the system escapes to the other local minimum on the left and so forth, and generates the whole free energy landscape [39]. but with the implicit solvent model. The atomic charges were derived from electrostatic potentials (ESP [47]) calculation using Gaussian09 at the MP2/6-31G** level based on the RESP (restrained ESP) protocol [48; 49; 50; 51]. The force field parameters for [EMIM]Cl were improved based on previous studies [41; 44], and it has been validated in our earlier work [52].