*2.2. Simulation Details*

The initial crystal structure of catechol (identifier CATCOL17) was obtained from CSD (Cambridge Structural Database), derived from experiments by Clegg and Scott, with the space group of P21/n (*a* = 9.8326 Å, *b* = 5.5910 Å, *c* = 10.467 Å, α = γ = 90◦, β = 114.988◦) [34]. As shown in Figure 1, there are four catechol molecules in the unit cell.

The entire geometry optimization and MD simulation process were implemented using the program Materials Studio 6.0 (Accelrys, San Diego, CA, USA) [35] with the COMPASS force field, which is a powerful ab initio force field with many accurate predictions of structures [36]. After the initial downloaded crystal was optimized using the COMPASS force field, the similar lattice parameters

of the two crystals listed in Table 1, proved that it was a suitable choice to use the COMPASS field since the relative error of each parameter was less than 5%.

**Figure 1.** The molecular structure (**a**) and the unit cell (**b**) of catechol.

**Table 1.** The comparison of the initial and the optimized lattice parameters of catechol.


After geometry optimization, the crystal habit of catechol in vacuum was predicted with the growth morphology method in the morphology module, providing several potential growing crystal faces that possibly existed in the solution environment. Based on the acquired morphologically important faces (h k l), the molecular models of catechol crystal surfaces were built by cleaving the corresponding (h k l) faces from the crystal cell with a depth of 2–4 *d*hkl and then extending the surface to supercells containing m × n unit cells (specific data for each crystal face listed in Supplementary Materials Table S1). The solvent effects were modeled using the Amorphous Cell module by setting up a three-dimensional solvent box containing 200 randomly distributed solvent molecules and the size of the solvent box was consistent with the crystal supercell and the thickness determined by the solvent density.

Subsequently, a double-layer crystal–solvent interface model was constructed, and then a vacuum layer with a thickness of 50 Å was added onto the model box to avoid any additional effects of periodic boundaries. After the interface model was geometrically optimized with catechol molecules fixed in the crystal, the MD simulation was operated with NVT ensemble at 298 K, controlled by the Andersen thermostat. The total duration of the simulation was 200 ps with a time step of 1 fs, and the data were collected every 100 steps. In the non-bonding interaction calculation procedure, the standard Ewald method was used to calculate the electrostatic interactions with an accuracy of 0.0001 kcal mol−<sup>1</sup> while the atom-based summation method was applied to calculate van der Waals interactions with a cut-off distance (*d*c) of 15.5 Å. The cleaving and extending operations were implemented to balance the accuracy and simulation time [37], assuring that the length and width of the crystal supercell were both no less than 2*d*c, while the thickness of the solvent layer and crystal layer were greater than *d*c [38].
