Ab initio determination of crystal stability of di-p-tolyl disulfide.

With the rapid growth of energy demand and the depletion of existing energy resources, the new materials with superior performances, low costs and environmental friendliness for energy production and storage are explored. Di-p-tolyl disulfide (p-Tol2S2) is a typical lubricating material, which has been applied in the field of energy storage. The conformational properties and phase transformations of p-Tol2S2 have been studied by pioneers, but their polymorphs and the polymorphism induced crystal structure changes require further analysis. In this study, we perform the crystal structural screening, prediction and optimization of p-Tol2S2 crystal with quantum mechanical calculations, i.e., density functional theory (DFT) and second-order Møller-Plesset perturbation (MP2) methods. A series of crystal structures with different molecular arrangements are generated based on the crystal structure screening. As compared to long-established lattice energy calculation, we take an advantage of using more accurate technique, which is Gibbs free energy calculation. It considers the effects of entropy and temperature to predict the crystal structures and energy landscape. By comparing the Gibbs free energies between predicted and experimental structures, we found that phase α is the most stable structure for p-Tol2S2 crystal at ambient temperature and standard atmospheric pressure. Furthermore, we provide an efficient method to discriminate different polymorphs that are otherwise difficult to be identified based on the Raman/IR spectra. The proposed work enable us to evaluate the quality of various crystal polymorphs rapidly.

Introduction

Organic class="Chemical">disulfide has been one of the major interesting research materials due to its suclass="Chemical">perior lubricating class="Chemical">proclass="Chemical">perty[1] and the extensive aclass="Chemical">pclass="Chemical">plications in storing energy[2]. class="Chemical">pan class="Chemical">Di-p-tolyl disulfide (p-Tol2S2, with the formula of C14H14S2), a typical representative molecule of a large group of diphenyl disulfides, has been explored in many aspects[3-5]. A number of analyses have examined the interplay of conformational transformation and the thermodynamic stability of the p-Tol2S2, demonstrating that the change in pressure[6-9] or recrystallization[10-12] would result in structural phase transition. Previous studies have primarily concentrated on the pressure effect on the conformational qualities of p-Tol2S2 and its molecular aggregation in crystalline state[13-17], but the polymorphs and polymorphism induced difference either in the crystal structure or the physical and chemical properties have not been studied. The present paper focuses on the polymorphs of p-Tol2S2, studying the stability of different forms and providing an effective spectral method to distinguish them. There are several different polymorphs of p-Tol2S2 with different conformations and molecular aggregations. The two main structures are phase α and phase β. Phase α is an independent molecule (Z′ = 1), which has an asymmetric conformation, with the monoclinic unit cell of space group P21. At higher pressures (> 1.6 GPa), phase α transform into a steric hindrance and aggregated phase β with reduced lattice parameters and volume. Structure β is a triclinic phase with two independent molecules (Z′ = 2) per unit cell and has a space group of P1 (see Fig. 1). The structures of p-Tol2S2 were determined by X-ray diffraction and its cell parameters can be found in the Cambridge Structural Database (CSD).

Figure 1

Crystal structures of p-Tol2S2. Phase α (a) has a space group of P21, while phase β (b) has a space group of P1. The molecules are colored by the symmetry operation, where the unit cells of phase α and phase β contain one and two independent molecules, respectively.

Crystal structures of class="Chemical">p-Tol2S2. class="Chemical">pan class="Chemical">Phase α (a) has a space group of P21, while phase β (b) has a space group of P1. The molecules are colored by the symmetry operation, where the unit cells of phase α and phase β contain one and two independent molecules, respectively.

class="Chemical">Polymorclass="Chemical">phisms, differing in their class="Chemical">physical class="Chemical">proclass="Chemical">perties, have different solubility and dissolution rates. Therefore, discrimination of class="Chemical">polymorclass="Chemical">phisms is imclass="Chemical">portant for all molecular crystal class="Chemical">products and industrial class="Chemical">production. Crystal structure class="Chemical">prediction (CSclass="Chemical">pan class="Chemical">P) is an emerging tool to help crystallographers determine various thermodynamically favorable crystal packings which are independent of the kinetics of crystallization. However they depend on a wide range of controllable or less easily controllable crystallization condition[18-20]. The advancement and development of CSP have been widely performed on different polymorphisms[21-23]. In this paper, to obtain more reliable structures, we use different CSP tools of MOLPAK[24,25] (MOLecular PAcKing) and USPEX[26-28] to perform crystal structure screening, prediction and optimization of p-Tol2S2. The MOLPAK and USPEX totally produce 9000+ possible hypothetical crystal structures which lie in the most common space groups with different molecular arrangements. The 38 candidates (20 from MOLPAK and 18 from USPEX) with the lowest lattice energies are selected to perform further optimization, where the crystal structure optimization is calculated at the DFT level using the ωB97XD/6-31G* method, and the single-point electronic energy is calculated at the MP2 level on the optimized crystal structure. The MP2 calculation is used to obtain more accurate and reliable prediction results. Besides, to treat the macromolecular system, such as the p-Tol2S2 (has 30 atoms) crystal in the present paper, the embedded-fragment method[29-31] is used to obtain the total energy by taking an appropriate combination of the energies of monomers and dimers, which are embedded in the electrostatic field of the rest in the crystalline environment. Via the Gibbs free energy calculation, we determine that phase α is the most stable structure from those predicted and experimental structures.

Methods

Crystal structure predictions

The credibility of computing methods for crystal structure prediction has been improved over the past years. Here we use the MOLclass="Chemical">PAK and USclass="Chemical">pan class="Chemical">PEX tools to determine the crystal forms from the molecular structure. MOLPAK can create unit cells that have the smallest volumes per molecule by allowing surrounding molecules to collapse along the three axes of a Cartesian coordinate system until a predetermined repulsion threshold with the central molecule is reached[32,33]. While the USPEX program, adapting the evolutionary algorithm that can find the stable structures for systems with various rigid molecular shapes[26,34]. The initial structures are generally generated by selected space groups with the initial known structures under atmospheric pressure. Structure relaxation is implemented from low to high precision using the VASP program. All the structures are examined and their initial sorting was performed based on lattice energy or enthalpy after full relaxation. However selected potential structures were resorted according to root-mean-square deviations (RMSDs), given in Tables S1 and S2.

Structure optimization and energy calculation

After generation of the initial crystal structures using MOLclass="Chemical">PAK and USclass="Chemical">pan class="Chemical">PEX, for further calculation, 20 structures from MOLPAK and 18 structures from USPEX were selected based on lower lattice energies and enthalpies, respectively. To minimize the enthalpy, we adopt quasi-Newton algorithm[35] for the crystal structure optimization. Because of the large molecular size and the quantity of molecules in the supercell, the embedded-fragment method[30,31,36] with ωB97XD/6-31G* was used for the calculation of enthalpy and Gibbs free energy of the p-Tol2S2 crystal. The Hessian matrix approximation was updated using the BFGS procedure[37] and the convergence criterion for the maximum gradient was set to 0.001 Hartree/Bohr. The embedded fragment method, as the name implies, is a process that breaks a large system into small fragments. The embedded fragment method divides the total energy per unit cell of the crystal into an appropriate combination of the energies of monomers and dimers, and therefore can treat the macromolecules effectively. The individual molecule is designated as a segment, and the interaction energies between two segments, in close contact with each other, are calculated by quantum mechanics (QM). The interaction energy between the two distal segments is calculated according to Coulomb interactions. The DFT calculations, accounting for the environmental influences, are embedded in an electrostatic field represented by the point charges of the remaining system.

The internal energy of a unit cell for the molecular crystal system is calculated by:where the three-integer index of unit cell is represented by a variable ‘n’. The class="Chemical">quantum mechanical (class="Chemical">pan class="Chemical">QM) energy of the ith molecules in the nth unit cell is represented by . The QM energy of dimers is represented by , where ith and jth represent number of molecules with respect to 0th central unit cell and nth unit cells[29,38-41]. The crystal system is represented by a supercell. The 1st part of Eq. (1) calculates all the single molecular energies in the 0th unit cell, which is in the center. The 2nd part represents the two-body QM interaction that has a shorter distance than (where is a given cutoff distance which is set to 4 Å in this work). The 3rd part of Eq. (1) provides the interactions between one molecule in the central unit cell and the other in nth unit cell whose distance is shorter than . The QM is used to calculate short range interactions (the first three parts in Eq. (1)). It is calculated in electrostatic field of the rest, where ωB97XD/6-31G* level is used to fit electrostatic potential variations. The background charges are represented by the supercell. The Coulomb’s charge-charge interaction is employed to approximately treat long-range interactions between two molecules of dimers whose distance is larger than . The long-range electrostatic interactions are represented by in a supercell.

The eclass="Chemical">nthalclass="Chemical">py for each unit cell can be calculated by considering the effect of external class="Chemical">pressure as follow:where the external class="Chemical">pressure is reclass="Chemical">presented by class="Chemical">pan class="Chemical">P and the unit cell volume is represented by V. The harmonic approximation is used to calculate zero-point vibrational energy and the entropy , which are shown in Eqs. (3) and (4), respectively, where the frequency of phonon with lattice vector k is represented by . and is the Boltzmann constant. The capital K is the product of all k, which are evenly spaced grid points in the reciprocal unit cell. The k-grid of 21 21 21 has been used in this study, where K = 9261. The calculation of Gibbs free energy with effects of temperature and pressure in a unit cell is given by Eq. (5).

As class="Gene">MP2/6-31G* is more accurate than DFT/ωB97XD/6-31G* but its comclass="Chemical">putational time and cost are very high. Therefore, to be more efficient and effective, we oclass="Chemical">ptimized crystal structure as well as calculated eclass="Chemical">pan class="Chemical">nthalpy (H) over single point energy, zero-point vibrational energy (), and entropy () at DFT/ωB97XD/6-31G* level. Moreover, to obtain better accuracy, MP2/6-31G* level was used to recalculate enthalpy (H) only, over single point energy based on crystal structure calculated by DFT/ ωB97XD/6-31G*.

Raman spectra simulation

For a periodic molecular system, the dynamical force constant matrix is calculated as follow:where the masses of atoms p and pan class="Chemical">q are reclass="Chemical">presented by m and m, resclass="Chemical">pectively.

The k is a Brillouin zone point which is known. The 2nd order derivative of the total energy in a unit cell which correspond to atom p and class="Chemical">q in the 0th and the class="Chemical">pan class="Chemical">nth cell, is represented by H(rp,0, rq,n) in Eq. (6) at geometry equilibrium. The number of nearby unit cells for QM treatment was truncated at R = 1, and the number of k-points was set to 21 in each of the three dimensions. The vibrational frequencies and the corresponding normal modes of the periodic molecular system can be obtained by solving the dynamical matrix . The zone-center (k = 0) vibrations presenting nonzero intensities results in the activation of IR- and Raman spectra. The force-constant matrix D(0) was used to calculate the vibrational frequency. Therefore, we can obtain the IR intensity (I) and Raman intensity (R) for the fundamental transition of the mode in the nth phonon branch with wave vector k = 0, demonstrated by the following derivatives:

The corresponding normal mode is represented by pan class="Chemical">Q. The diclass="Chemical">pole moment derivative and the class="Chemical">polarizability derivative in the central unit cell can be derived with the embedded-fragment class="Chemical">pan class="Chemical">quantum mechanical method.

Results and discussion

The experimental data of the class="Chemical">p-Tol2S2 (CH3-C6H4-S)2 in the Cambridge Structural Database[42,43] is used to class="Chemical">predict the crystal structure. class="Chemical">pan class="Chemical">Phase α has a space group of P21, with a monoclinic unit cell and a volume of 628.741 Å3. The cell lattice parameters are a = 7.5927 Å, b = 5.71318 Å, c = 14.722 Å, α = 90°, β = 94.7615°, γ = 90°. However, phase β has a space group of P1 with a unit cell volume of 545.963 Å3. The lattice parameters of unit cell are a = 7.3057 Å, b = 5.5093 Å, c = 14.038 Å, α = 95.14°, β = 97.23°, γ = 85.36°. Figure 1 represents the molecular structures of both (α and β) phases.

The crystal structure prediction was performed by the MOLclass="Chemical">PAK and USclass="Chemical">pan class="Chemical">PEX programs. From MOLPAK, we obtained 31 standard clusters of crystal structures with common space groups of (P212121, P, P21/c, Cc, C2/c, P21, Pna21 and Pbca). Each cluster contains 300+ crystal structures. Therefore, MOLPAK generated total crystal structures more than 9000. From these 9000+ crystal structures, only 20 crystal structures were selected according to lower lattice energies, and these selected structures were further optimized. All the selected 20 structures from MOLPAK have same lattice energy of − 16.91 kcal/mol, given in Table S1. We obtained same lattice energy for 20 selected structures because of two main reasons: first, these structures are similar to each other with minor difference in their lattice parameters, and second, the MOLPAK software does not calculate lattice energy very accurately. Therefore, to obtain accurate energy landscape, we optimized selected structures and calculated their Gibbs free energy at the QM level.

On the other hand, USclass="Chemical">PEX generated 20 generations with each generation generated about 20 crystal structures. Therefore, the total crystal structures generated by USclass="Chemical">pan class="Chemical">PEX were about 400. From USPEX predictions, only 18 crystal structures were selected based on enthalpy sorting. After selection of potential candidates we sorted them again, according to RMSDs, given in Table S2. The lattice parameters of phase α, phase β and the 38 predicted new crystal structures by MOLPAK and USPEX are given in Tables S1 and S2 as Supporting Information.

The Gibbs free energy

The stability of 40 crystal structures are compared in Fig. 2, including phase α, phase β and 38 predicted new structures by MOLclass="Chemical">PAK and USclass="Chemical">pan class="Chemical">PEX. Figure 2 shows the difference of Gibbs free energy per unit cell among 38 predicted structures, phase α and phase β of the p-Tol2S2 as a function of unit cell RMSD. In Fig. 2, the purple and green circles represent the 38 predicted structures by MOLPAK and USPEX, whereas the green square and red diamond denote the experimental Gibbs free energies of phases α and β, respectively. The enthalpy over single point energy of 40 crystal structures was calculated at the MP2/6-31G* level based on the DFT/ωB97XD/6-31G* optimized crystal structures.

Figure 2

The Gibbs free energy differences per unit cell between the predicted crystal structures and experimental phase α as a function of the unit cell RMSD (with respect to phase α) under 298 K and standard atmospheric pressure, where the purple and green circles signify the 20 MOLPAK and 18 USPEX predicted structures, respectively. The green square and red diamond represent Gibbs free energies of the phases α and β, respectively, which were obtained from experiment.

The Gibbs free energy differences per unit cell between the predicted crystal structures and experimental phase α as a function of the unit cell RMSD (with respect to phase α) under 298 K and standard atmospheric pressure, where the purple and green circles signify the 20 MOLclass="Chemical">PAK and 18 USclass="Chemical">pan class="Chemical">PEX predicted structures, respectively. The green square and red diamond represent Gibbs free energies of the phases α and β, respectively, which were obtained from experiment.

The values of RMSD vary from 1 to 3.5 Å, which suggest that the crystal structure prediction method is reliable. Moreover, the lowest value of RMSD is smaller than 2 Å demonstrating the accuracy of prediction. The correlation between RMSD and Gibbs free energy (as shown in Fig. 2) reveals that the free energy of phase α, locating at the bottom of the lattice energy landscape, is the most stable structure among all predicted and experimental structures. Furthermore, the selected 20 structures from MOLpan class="Chemical">PAK were class="Chemical">possessing same lattice energies, given in Table S1. However, in Fig. 2 Gibbs free energy varies from about 1.5 to 8 kcal/mol, which is evidence of accuracy of Gibbs free energy calculation.

To verify the rationality of predicted structures, we selected eight candidates from 38 predicted structures by USclass="Chemical">PEX and MOLclass="Chemical">pan class="Chemical">PAK programs to visualize overlay only. Figures S1 and S2 in the Supporting Information present the similarity of crystal structures and the molecular conformations via the calculations of the structural RMSDs, phase α and the predicted structures based on the RDKIT tool[44]. We selected a subset of 8 structures from 38 structures based on difference in RMSDs. The green molecules represent the structures of p-Tol2S2 phase α determined by X-ray diffraction and the red molecules denote the predicted structures based on USPEX and MOLPAK programs. The molecular overlay shows the reliability of the predicted structures.

Vibrational spectra

The crystal structure of a molecule can be identified by vibrational spectroscopy[45,46]. The calculated full range of Raman spectra are provided in Figs. S3 and S4 of the Supporting Information. For the purpose of highlighting the differences of phase α and phase β, we select parts of the calculated Raman and IR spectra of class="Chemical">p-Tol2S2 crystals. The Raman and IR sclass="Chemical">pectra are shown in Figs. 3 and 4. In Fig. 3, the Raman sclass="Chemical">pectra are calculated from 700 to 1200 cm−1, where class="Chemical">phase α (green curve) and class="Chemical">phase β (red curve) are class="Chemical">predicted to disclass="Chemical">play seven and eight discernible intense Raman bands, resclass="Chemical">pectively. The band 3*, locating at ~ 874 cm–1, in the red curve is the characteristic class="Chemical">peak for class="Chemical">phase β. Similarly, there are five and four IR bands for class="Chemical">phase α (green curve) and class="Chemical">phase β (red curve), resclass="Chemical">pectively. The band 2*, locating at ~ 3220 cm−1, in the green curve is the characteristic class="Chemical">peak for class="Chemical">phase α. Therefore, the band 3* in the Raman sclass="Chemical">pectra and the band 2* in the IR sclass="Chemical">pectra, are the characteristic bands to distinguish class="Chemical">phases α and β of class="Chemical">pan class="Chemical">p-Tol2S2. Figures 3 and 4 show that although the structures of phases α and β are very similar (as shown in Fig. 1), they can be identified by the characteristic bands in the vibrational spectra. The vibrational spectra of different pressure about phase α and phase β are demonstrated in Figs. S5–S8 in the Supporting Information. We can observe the vibrational spectra as a function of different pressures.

Figure 3

Calculated Raman spectra of p-Tol2S2 phase α (green) and phase β (red) under standard atmospheric pressure.

Figure 4

Calculated IR spectra of p-Tol2S2 phase α (green) and phase β (red) under standard atmospheric pressure.

Calculated Raman spectra of pan class="Chemical">p-Tol2S2 class="Chemical">phase α (green) and class="Chemical">phase β (red) under standard atmosclass="Chemical">pheric class="Chemical">pressure.

Calculated IR spectra of pan class="Chemical">p-Tol2S2 class="Chemical">phase α (green) and class="Chemical">phase β (red) under standard atmosclass="Chemical">pheric class="Chemical">pressure.

Conclusions

It is highly desirable to control polymorph for all the manufacturing industries which manufacture and/or utilize crystalline products. This will enable us to enhance functionality of the products to next higher level. A high precision method to screen different polymorphs for designing of crystals and to find out the most stable structure that confers an optimal product would, therefore, be extremely valuable. In this paper, we carried out a theoretical approach to screen different polymorphs of class="Chemical">p-Tol2S2, a lubricating and energy storage material, based on the DFT and class="Chemical">pan class="Gene">MP2 methods with the embedded fragment approach. With the theoretical calculation, we pick out the most stable structure phase α from the predicted and experimental structures. Considering the small RMSDs, there will be hundreds of predicted structures between them and the conformation of existing phases α and β, the proposed method has successfully performed the predictions for sufficiently similar conformations along with the identification of the most stable structure. The predicted vibrational spectra provide an effective way to identify polymorphs with the characteristic peaks of IR and Raman spectra.

Supplementary Information.