**1. Introduction**

Crystal size and shape have essential effects on downstream processing such as filtration, washing and drying for solution crystallization [1–3]. Meanwhile, crystal properties such as flowability and bulk density [4,5] are closely associated with crystal morphology. Therefore, it has been very crucial to explore the possible crystal habits under different conditions to select the suitable crystals for industrial operations. Many factors contribute to affect crystal shapes, such as solvents [6], additives [7], temperature [8], supersaturation [9] and even stirring rate [10], which can be modified and controlled to obtain the morphology optimal for industrial applications. For solution crystallization, the effect of solvents is one of the most primary factors that affects final crystal habits. As indicated by many researches [11–13], the interactions between solvent molecules and crystal faces are quite essential to crystal growth. Hence, in order to obtain desired crystal shapes, it is necessary to investigate solvent effects on crystal habits and aspect ratio changes.

Although real and effective crystal morphology could be acquired by experimental methods, it is labor-intensive and time-consuming to conduct a large number of experiments just to find the optimal one that benefits production. In recent years, computer-aided methods such as molecular simulations and first principles (FP) simulations have been explored to explain and predict crystal habits [14–17], and the gradually developing simulation methods have provided a broader perspective for crystal morphology research at a molecular level.

Molecular dynamics (MD) is one of the molecular simulation methods that helps to investigate the directions of crystal growth with molecular information, reveal the interactions between solvent molecules and crystal surfaces and provide more microscopic details for experiments. The MD method has been successfully applied to simulate crystal morphology in many cases [18–20]. Wang et al. [21] utilized a surface docking model for a crystal growth simulation of cefaclor dihydrate and proposed a competitive relationship of surface adsorption by the solute and solvent molecules based on pure solvent and solution models. Zhang et al. [22] simulated the crystal morphology of ibuprofen obtained from four different solvents using the modified attachment energy (MAE) model and found that the crystal aspect ratios were sensitive to the relative polarity of the solvents. They regarded the method as a promising way for the computer-aided design of desirable pharmaceutical crystal habits with rapid solvent screening. Poornachary et al. [23] investigated the mechanism of additive inhibition on naproxen crystal growth by modeling the intermolecular interactions between the polymeric additive and the crystal surface and ascribed the phenomenon to solute diffusion control.

However, the exploration of solvent–crystal interactions from the thermodynamic perspective is just a partial picture, and discussions are not sufficient on the differences in crystal aspect ratios because of solvent system distinctions, which makes the crystal morphology selection in industrial application still based on experiences.

Catechol is a common fine chemical raw material, which is mainly used as an intermediate for polymerization inhibitor synthesis and pesticide production. In addition, it is also consumed as a precursor to fragrances, pharmaceuticals and dyes [24]. This organic chemical commodity has been in industrial production for over forty years, and the research emphasis is usually focused on its synthesis techniques [25,26]. Further experimental and simulation investigations on catechol crystal morphology are needed to get a deeper insight into the crystal growth and habits in solvents, which helps to select optimal crystal shapes and solvents to avoid agglomeration and improve product quality for industrial production.

In this paper, the catechol crystal morphology in isopropanol, methyl acetate and ethyl acetate were simulated using the MAE model, which helped to quantify the interactions between the crystal surfaces and solvent molecules. Experiments of cooling crystallization were conducted at the same time for model verification. This study aimed to describe the solvent effects on crystal habits from the perspective of crystal–solvent interactions and surface properties, which may favor a better understanding of crystal habit distinctions, especially the aspect ratio changes of crystal shapes in solvents.

#### **2. Calculation Methodology**

#### *2.1. Theory*

The Bravais–Friedel–Donnay–Harker (BFDH) model [27] is one of the models which are initially applied to predict crystal habits, but the model lacks precision because it simulates possible crystal facets merely according to geometric factors without considering the actual chemical environment, so it was soon developed into the attachment energy (AE) model by Hartman and Bennema [28–30], taking into account the energies of the system on the basis of the period bond chain (PBC) theory [31]. This kind of model is based on the attachment energy (*E*att), which is defined as the released energy

when a growth slice attaches on to a growing crystal surface. The relative growth rate (*R*hkl) of each crystal face (h k l) is assumed to be proportional to the absolute value of the attachment energy [28]:

$$|R\_{\rm hkl}| \propto |E\_{\rm att}|.\tag{1}$$

The attachment energy (*E*att) is calculated as

$$E\_{\text{att}} = E\_{\text{latt}} - E\_{\text{slice}} \tag{2}$$

where *E*latt is the lattice energy of the crystal and *E*slice is the energy of a growth slice with a thickness of *d*hkl.

However, the AE model was put forward to simulate ideal crystal morphology in the premise of a vacuum environment, leading to a precision loss when the model was applied to solution crystallization. As a matter of fact, solvent molecules adsorb on the crystal surfaces as solute molecules do [32], which means that the growth of the crystal faces may be impeded by the solvent molecule adsorption. In order to cover the effects of the solvent on crystal growth, a modified attachment energy (MAE) model was developed. The modified model introduced an energy correction term *E*s for the initial *E*att to take into consideration the effects of the solvent molecules on crystal faces, and the modified attachment energy *E'*att is calculated in the following formula:

$$E'\_{\rm att} = E\_{\rm att} - E\_{\rm s.} \tag{3}$$

where *E*<sup>s</sup> represents the binding energy between the solvent molecules and the crystal face (h k l), which it can be obtained using the follow equation:

$$E\_{\rm s} = E\_{\rm int} \times \frac{A\_{\rm acc}}{A\_{\rm box}},\tag{4}$$

where *A*box is the total crystal surface area of the simulated supercell along the (h k l) direction, and *A*acc represents the solvent-accessible area of the crystal face in the unit cell, which can be approximated by the Connolly surface algorithm [33]. *E*int was defined as the interaction energy between a specific crystal face and the corresponding solvent layer, which can be calculated as follows:

$$E\_{\rm int} = E\_{\rm tot} - (E\_{\rm sur} + E\_{\rm sol}).\tag{5}$$

In Equation (5), *E*tot means the total energy of the entire simulation box including all crystal and solvent molecules, and *E*sur and *E*sol represent the energy of the isolated crystal surface and solvent layer, respectively.

After correction, the relative growth rate *R'*hkl in the MAE model was still in proportion to the absolute value of the modified attachment energy |*E'*att|:

$$|R'\_{\text{hkl}} \ll |E'\_{\text{att}}|.\tag{6}$$
