Identifying the key steps determining the selectivity of toluene methylation with methanol over HZSM-5

Methylation of toluene with methanol to produce p-xylene has been investigated for decades, but the origin of selectivity is still under debate. Here we report computational studies based on ab initio molecular dynamics simulations and free energy sampling methods to identify the key steps determining the selectivity. The steps of toluene methylation to protonated-xylene, deprotonation of protonated-xylenes, and diffusion of xylene in HZSM-5 channels are compared. We find the pathways of formation for protonated p-/m-xylenes have similar free energy barriers. Meanwhile, the methylation is found rate-determining, thus the probability to generate p-/m-xylenes at the active site are similar. We then find that the diffusion for m-xylene along the zigzag channel is more difficult than its isomerization to p-xylene, which in turn further promotes the selectivity of p-xylene formation. These insights obtained at the molecular level are crucial for further development of high-performance zeolite catalysts for toluene methylation.


Supplementary Note 1. Stepwise mechanism results
According to the stepwise mechanism as shown in Supplementary Figure 1, the methylation reaction was divided into two parts, i.e. adsorbed methanol dissociation to form surface methyl groups and surface methyl groups reacting with toluene to form proton-xylene.

Supplementary Figure 1.
Stepwise mechanism of toluene methylation reaction.
MTD method was used to study these reactions. As for the surface methyl formation reaction, CV1 was the CN of Ome-Cme, CV2 was the CN of Cme-Oz, nn, nd, r0 were 5, 10, 1.5 Å respectively (according to Equation 1 in the main text). The initial height of Gaussian hill was 5.251 kJ mol -1 . All other parameters are the same with the settings in the main text (details shown in Methods section). Supplementary Figure 2 shows that the free energy barrier is greater than 150 kJ mol -1 , which is close to the value of 152 kJ mol -1 reported in the literature. 1 One can find that the free energy barrier of surface methyl group formation is much higher than the effective barrier of the concerted mechanism, and therefore the stepwise mechanism should be preferred under the conditions considered in our work. Figure 2. The curve showing free energy barrier and reaction free energy of surface methyl groups formation against time.

Supplementary
As for the proton-PX/MX formation reaction, CV1 was the CN of Cme-Oz, nn, nd, r0 were 6, 12 and 1.5 Å, respectively; CV2 was the CN of Cme-Cp/m, nn, nd, r0 were 6, 12 and 2.0 Å, respectively. Other parameters remain unchanged. Supplementary Figure   3 shows the free energy barriers are quite similar and proton-PX is more stable than proton-MX in HZSM-5, which are consistent with the results of direct methylation mechanism. Figure 3. The curve of free energy barrier and reaction free energy of (a) proton-PX and (b) proton-MX formation reactions against time.

Supplementary Note 2. NPT optimization of lattice parameters
A 50 ps NPT simulation was performed at the temperature of 670 K and the pressure of 1 bar, to confirm the reliability of lattice parameters used in the current work.
By averaging over the range from 5 to 50 ps (the first 5 ps was used to equilibrate the system), the optimized lattice constants were evaluated to be a = 20.224 Å, b = 20.014 Å, c = 13.485 Å, which is very close to the values used in the original manuscript with a negligible discrepancy of ~ 0.02 Å (~ 0.1 %). and FS (proton-PX) states were obtained in the 2D free energy surface using a certain energy range. The energy range for the IS and FS region is from AIS/FS to AIS/FS + kBT, while for TS this is from ATS -0.5 kBT to ATS + 0.5 kBT (1 kBT = 5 kJ mol -1 at 670 K).

Supplementary
Then rectangle regions that meet the criteria were selected, as shown in Supplementary  In addition, Supplementary Figure 9 To confirm that the Oz-Hp/m distance of 3.0 Å can represent the free energy minimum, we set the reaction coordinate change rate to 0.000109 Å frame -1 , while the initial structure and other parameters remain unchanged. As for proton-PX, the free energy change is small and fluctuates around 0, and for proton-MX the free energy increases gradually, demonstrating that initial reaction coordinate approaches to the free energy minimum.

Supplementary Note 9. Diffusion barriers
The values of diffusion barriers and its uncertainties were shown in Supplementary   Table 5. The uncertainties of free energy barriers were obtained in the following way: Firstly, the standard deviation (sf) of constrained force was estimated using the block average method. Then the linear error propagation theory was used to calculate the uncertainty of free energy profiles by summing up the variance of work in each bin (sw 2 ) from minimum to maximum of the free energy profiles. Two examples trajectories of Fc were given in Supplementary Figure 12. From the trajectories of constraint force Fc, we can see that Fc is equilibrated in the range of 1000 to 5000 steps. Therefore, in each constrained MD simulation, totally 5000 steps were performed and the last 4000 steps of the trajectory were analyzed.