Vibrational Spectra and Molecular Vibrational Behaviors of Dibenzyl Disulfide, Dibenzyl Sulphide and Bibenzyl

The vibration spectroscopy (Raman and infrared) of widely concerned molecules in sulfur corrosion phenomenon (Dibenzyl Disulfide, Dibenzyl Sulphide, and Bibenzyl) is detailedly analyzed based on density functional theory and experimental measurement. The dominant conformations of these molecules are determined according to Boltzmann distribution in relative Gibbs free energy. Additionally, noncovalent interaction analysis is conducted to indicate intramolecular interaction. Vibration normal mode is assigned based on potential energy distribution, which comprehensively reveals the molecular vibrational behaviors. Conformations weighted spectra are obtained and compared with experimentally measured spectra. We found that experimental spectra are in good agreement with the theoretical spectra in B3LYP-D3(BJ)/6-311G** level with a frequency correction factor. Furthermore, the divergence among these molecules is discussed. The vibrational behavior of the methylene group in the molecule shows a trend with the presence of the sulfur atom.


Introduction
Nowadays, the sulfur corrosion problem of power equipment has received widespread concern [1][2][3][4][5]. The failures caused by corrosive sulfur greatly affect the safety and stability of the power system and result in substantial economic losses. Dibenzyl Disulfide (DBDS) is regarded as the most important corrosive sulfur. The corrosion mechanism was proposed by Mitsubishi Electric Corporation as shown in Figure 1 [6]. The intermediate product Dibenzyl Sulphide (DBS) is also corrosive, and the final product of this reaction is Bibenzyl (BiBz) and cuprous sulfide. The DBDS, DBS, and BiBz are considered as indicators in power equipment sulfur corrosion phenomenon. The study on their physical and chemical properties will contribute to the development of advanced corrosive resistant materials and the realization of corrosive state monitoring. Numerous efforts have been made to study properties of the DBDS, DBS, and BiBz, including chemical reactivity [7,8], synergistic effect [9,10], detection methods [6,11], etc. [12,13]. In this paper, the vibrational spectra and molecular vibrational behaviors are analyzed detailedly. The vibrational spectroscopy includes infrared and Raman spectroscopy, which is based on periodic changes of dipolmoments (IR) or polarizabilities (Raman) caused by molecular vibrations behaviors. IR and Raman spectra are both widely used methods for composition analysis and trace detection [14][15][16].
In addition to experimentally measured vibrational spectra, density functional theory (DFT) provides a solid foundation of theoretically calculated spectra. Considering the presence of a large number of fundamentals and overtones bands, it is difficult to interpret experimental spectra of large molecule [17]. However, the conformational distribution and theoretical spectrum of molecules can be obtained based on DFT. Then, a potential energy distribution (PED) based analysis can be performed to achieve vibration assignment, which gives a clear and crispy concept about the structural-vibrational properties of the molecule. A lot of studies have recently appeared in the literature concerning the calculation of spectroscopic properties by DFT methods [17][18][19][20][21][22].
A confocal Raman spectrometer mentioned in our previous paper [23] was used to detect experimental spectra of the DBDS, DBS, and BiBz. There are three kinds of gratings (600 lines/mm, 1200 lines/mm, and 1800 lines/mm) built in the spectrometer. The grating with minor grooves number is designed for wide spectral range detection, although the spectral resolution will decrease. In order to compare with theoretical spectra in full spectral range, 600 lines/mm grating was adopted. The other test parameters of the measurement are listed below: the slit width was 50 um, the integration time was 1 s, and the integration number was 20 times. The Raman platform was calibrated with monocrystalline silicon before each test. The Raman shift of the first-order silicon peak was adjusted to 520.7 cm −1 . The result of Raman spectrum comes from the average of three measurements.
The Raman/IR spectra of DBDS, DBS, and BiBz were plotted in Section 4 after normalization and de-baseline operations.

Conformation Determination
Flexible molecules always have multiple conformations, such as DBDS, DBS, and BiBz studied in this paper. The proportion of different conformations was calculated according to Boltzmann distribution in relative Gibbs free energy. The molecular conformation search program Molclus [24] is used for conformation determination in this paper.
Firstly, a batch of initial conformations of DBDS, DBS, and BiBz were generated by rotating chemical bond via Molclus. The chemical bond within the benzene ring was ignored; only C-C bonds, C-S bonds, and S-S bonds between benzene rings are considered. All the considered bonds rotate every 120°. The DBDS with five rotation bonds generates a total of 3 5 initial conformations. Similarly, DBS generates 3 4 initial conformations. BiBz generates 3 3 initial conformations.
All initial conformations generated above were pre-optimized geometrically by Gaussian [25] in a semi-empirical method PM7 [26]. The optimized conformations with the energy difference within 0.25 kcal/mol and the geometric difference within 0.1 Å will be classified as the same conformation. The conformations with relatively lower energy were applied further geometric optimization and frequency analysis in B3LYP-D3(BJ)/6-311G** [27][28][29][30][31] level. For an authentic conformation distribution, a double hybrid functional revDSD-PBEP86-D3(BJ) [32] in conjunction with may-cc-PVTZ [33] was introduced to recalculate the single point energy. The single point can be added with thermal correction to find the precise Gibbs free energy. The thermal correction is obtained in frequency analysis based on rigid-rotor harmonic oscillator (RRHO) model [34] under 298.15K and 1 atm. The results are shown in Table 1. Conformations in Table 1 with a proportion less than 5% were ignored. The geometric structures and surface electrostatic potential of several dominant conformations were visualized via VMD [35], as shown in Figure 2. The electrostatic potential involved was evaluated by Multiwfn, based on the highly effective algorithm proposed in ref. [36]. It can be found from the surface electrostatic potential distribution that the negative charge region mainly occurs around the sulfur atom and the benzene ring planar center. It is difficult to explain the geometric structure only by the surface electrostatic potential distribution due to the existence of various intramolecular interaction. Further NCI analysis is thus needed.

Noncovalent Interaction Analysis of the Dominant Conformation
To reveal the nature of these dominant conformations with low free energy, the noncovalent interaction (NCI) [37] method was applied for weak molecular interactions analysis.
The NCI method uses Reduced Density Gradient (RDG) to distinguish regions with noncovalent interactions in the system. The expression of RDG is as follows, where ρ is the electron density.
The noncovalent interactions areas can be determined by RDG value. However, the strength and type of interaction are determined by ρ and sign(λ 2 ), respectively. The sign(λ 2 ) is the second largest eigenvalue of the electron density Hessian matrix. Overlapping the color map of sign(λ 2 ) · ρ, as shown in Figure 3, on the RDG value isosurface, means the areas, intensity, and type of the noncovalent interactions can be visualized at the same time. The isosurface map and scatter plot of five dominant conformations were visualized via Multiwfn [38,39] and VMD [35]. As described in Figure 3f, the red area in the isosurface map and the scatter plot indicates the repulsion between molecules, such as steric effect. The green area indicates the van der Waals interaction, while the blue area indicates hydrogen bonds and halogen bonds that have a significant attraction.
There are two common features in the five dominant conformations. Firstly, significant steric effect is observed in the two benzene rings. This steric effect is shown as a red spindleshaped area in the benzene ring on the RDG isosurface map, and as a red spike on the right side in the scatter plot. Secondly, there is no notable attraction in these conformations.
For BiBz-1, the van der Waals interaction is formed between the two benzene rings. The molecular structure of BiBz-2 is stretched, the two benzene rings are far apart, and no significant van der Waals effect is observed. This was shown as two spikes near the zero point on the scatter plot. The spikes can be positive or negative as the sign of sign(λ 2 ) · ρ is relatively unstable as the electron density in the van der Waals area is small. The BiBz-1 formed by two benzene rings attracting each other is more stable and has a lower single point energy due to intermolecular forces. However, as BiBz-2 has a lower thermal correction, the Gibbs free energy and proportion of these two conformations are not remarkably different.
For DBS, in addition to the zone between benzene rings, the van der Waals effect also occurs between the sulfur atom and the hydrogen atom on the benzene ring. However, this kind of H-S interaction is distinguished from hydrogen bonds (shown in blue). The RDG isosurface between the two benzene rings in DBDS is much more complicated, but in conclusion, it is still van der Waals interaction. These broken isosurfaces were shown as many fine spikes near the zero point on the scatter plot in Figure 3d,e.

PED Based Vibrational Assignment
Each characteristic group in molecules has its own unique vibration behavior. However, these groups will be coupled with the molecular systems, which change frequency and intensity. Non-linear molecules have 3n−6 normal vibration modes, where n is the number of atoms in the molecule. The normal modes are the coupling of multiple groups vibration. It is difficult to directly identify which groups and vibration types contribute to the normal mode.
PED (potential energy distribution) analysis is an effective method that decomposes the normal vibration mode into group characteristic vibrations. The contribution proportion of the group vibration can be obtained, and the vibrational features of normal mode can be revealed. In this paper, PED based vibrational assignment of DBDS, DBS, and BiBz was performed via VEDA4 program [40].
The non-linear molecule DBDS consists of 30 atoms, it has a total of 84 normal modes. Therefore, 84 linearly independent internal coordinates are introduced to describe molecular vibration. The internal coordinates and related atoms of the DBDS-1 conformation are shown in Table 2. Atom numbers are listed in Figure 4. Additionally, the internal coordinates and related information of other four conformations can be found in Supplementary Material. The vibration type in VEDA4 program are divided into stretching (ν), bending (δ), and torsion (τ/γ), which correspond to specific internal coordinates. The stretching vibration involves two atoms. Bending vibration involves three atoms and torsion vibration involves four atoms. Among them, torsion can be further divided into two types: τ and γ. In many cases, the two definitions can describe the same atoms movement, but not always. τ ABCD means the dihedral angle between the ABC and BCD planes. γ ABCD means the angle between the AD vector and the BCD plane (reffed from VEDA4 user document). Table 2. Internal coordinates and related atoms in DBDS-1.

Stretching
Related Atoms Bending Related Atoms Torsion Related Atoms  The internal coordinates of DBDS-1 conformation consist of 29 kinds of stretching, 28 kinds of bending, and 27 kinds of torsion. The contribution of different internal coordinates to the normal vibration mode is shown in Tables 3 and 4. Similarly, the detailed table of other conformations is provided in Supplementary Material. It should be noted that the internal coordinates contribution less than 10% are ignored. The normal modes of DBDS-1 are divided into low frequency (<400 cm −1 ), middle frequency (400∼2800 cm −1 ), and high frequency (>2800 cm −1 ). The scaled frequency was corrected by a factor of 0.9640 [41].   Modes #1∼#14 are high frequency normal modes. All of them are contributed by stretching vibration, and there is a large frequency gap between #14 and #15. They have completely different characteristics in the frequency band. This is because the force constant of the stretching vibration is larger. The stretching vibration energy is greater than bending and torsion vibration energy. Therefore, the characteristic frequencies corresponding to the stretching vibration are all located in the high frequency region. Generally, almost all molecules have abundant vibration peaks corresponding to stretching vibrations in the high frequency region. The main feature of the low frequency zone is that no stretching contribution is assigned, while the middle frequency area is composed of various vibration contributions. The spectrum characteristics in the low frequency region is ignored in the following analysis. As the measurement range of FTIR only covers 400-4000 cm −1 , and the range of Raman spectroscopy in the low frequency zone is seriously affected by Rayleigh scattering.
The normal modes always appear in pairs. This is because DBDS is a quasi-symmetric molecule. Although the DBDS-1 conformation does not have molecular symmetry, its molecular composition is still symmetrical. For example, the vibration mode of one benzene ring must have the corresponding vibration mode of another benzene ring.

Methylene(-CH 2 ) Group
In the high frequency zone, the asymmetric stretching vibration of the two -CH2 groups corresponds to the #11 and #12 normal modes. The two mode peaks appear at 3007 cm −1 and 3025 cm −1 respectively, due to the vibration coupling (the calculated Raman and IR spectra are shown in Figures S4 and S5, respectively). The infrared activity of #12 mode is weak, almost invisible, in the IR spectrum. The symmetric stretching vibration of the two -CH2 groups corresponds to the #13 and #14 normal modes. These two modes appear as the same peak at 2946 cm −1 . The frequency of -CH2 symmetrical vibration is often smaller than that of asymmetrical vibration [21].
In the middle frequency zone, the scissoring vibration of -CH2 appears at 1416 cm −1 , corresponding to mode #23 and #24. When wagging vibration of -CH2 appears at 1223 cm −1 , corresponding to mode #29 and #30. High Raman activity and infrared activity are simultaneously exhibited in wagging vibration. The twisting vibration of -CH2 corresponds to #37 and #38 at 1126 cm −1 . The rocking vibration of -CH2 corresponds to #51 and #52 at 863 cm −1 . The activity of these four modes is relatively weak.
For the -CH2 group, the order of normal vibration frequency from large to small is: asymmetric stretch, symmetric stretch, scissoring, wagging, twisting, and rocking.

Phenyl(-C 5 H 6 ) Group
In the high frequency zone, the #1 and #2 normal modes correspond to the in-phase stretching vibration of the C-H bond on the benzene ring. For both #1 and #2 modes, the five C-H bonds on the benzene ring will vibrate in the synchronous phase. Two modes appear as the same peak at 3077 cm −1 . This characteristic peak is the strongest peak in the Raman spectrum, but not in the infrared spectrum. The eight #3∼#10 normal modes correspond to the asynchronous phase vibration of the C-H bonds on the benzene ring, and the Raman activity and infrared activity of these eight modes are quite different. Modes #3 and #4 appear at 3069 cm −1 in the IR spectrum. Modes #5 and #6 appear at 3062 cm −1 in both the IR and Raman spectrum. Modes #7 and #8 appear as a peak at 3051 cm −1 in the Raman spectrum.
In the middle frequency zone, the stretching vibration of the carbon skeleton of the benzene ring corresponds to modes #15∼#18, #27, and #28. The frequency of these stretching vibration modes is much smaller than the C-H bond stretching vibration. This shows that the vibration force constant of the carbon skeleton is much smaller than the force constant of the C-H bond. The bending vibration of the hydrogen atoms on the benzene ring corresponds to #19∼#22, #25, #26, #33∼#36, #39, and #40. Modes #31 and #32 have very strong coupling characteristics. All atoms except disulfide bonds are involved, and both apparent carbon skeleton stretching and hydrogen atoms bending occurred.
The #41 mode (frequency < 1016 cm −1 ) began to exhibit bending of the benzene ring carbon skeleton The bending and stretching vibrations of the carbon skeleton of the benzene ring simultaneously appear in the #41 and #42 modes. The #43, #44, #55, #56, and #63∼#66 modes correspond to bending vibration of benzene ring carbon skeleton. High Raman activity is exhibited in these two modes, appearing as peak at 986 cm −1 , which is also the characteristic peak used for Raman spectroscopy based DBDS detection in our study in progress.
From the #45 mode (frequency < 967 cm −1 ), the out-of-plane bending of hydrogen atoms on benzene ring occurs. Modes #45∼#50, #53 and #54, and #57∼#60 correspond to this type of vibration. Modes #67∼#69 correspond to the deformation vibration of the benzene ring couple with the stretching vibration of disulfide bond. The other out-of-plane vibration of the benzene ring appeared in the low frequency zone.
The #61 and #62 modes are composed of multiple vibrations, mainly including the in-plane bending, out-of-plane bending of the benzene ring, and the stretching vibration of C-S bond.

Weighted Vibrational Spectra
Different conformations have their own vibrational spectra. The weighed spectrum was calculated by different weighted conformations spectra according to their proportion. Frequency calculation is carried out in the optimal conformations; there is no imaginary frequency observed in the result. A frequency correction factor (0.9640) for B3LYP/6-311G** [41] was applied to the vibration spectra of DBDS, DBS, and BiBz due to the overestimated frequency of anharmonic effect.
Take the Raman spectra of DBDS as an example to illustrate the role of weighted spectrum. Figure 5 shows the Raman spectra of two DBDS conformations in the high frequency zone. The most significant difference between the two conformations is that DBDS-1 has a characteristic peak at 3025 cm −1 , but DBDS-2 has a characteristic peak at 2970 cm −1 .
Based on the vibration assignment in Section 4.1, the peak at 3025 cm −1 is assigned to asymmetrical stretching vibration of two -CH 2 groups in DBDS-1. The peak at 2970 cm −1 is assigned to symmetric stretching of -CH 2 group in DBDS-2. The Raman spectra of different conformations are different, but they are all superimposed in the weighted curve, as shown in Figure 5's blue line. The vibration spectroscopy is revealed much better when based on conformation weighted curve. It should be noted that the Raman spectrum obtained by theoretical calculation is a spectrum of Raman activity (unit in Å 4 /amu). In order to compare with the experimental spectrum, the Raman activity is converted into the Raman intensity, depending on the frequency of incident light and ambient temperature [42]: where S i , I i , and v i are the Raman activity, Raman intensity, and vibration frequency of ith vibration mode respectively. v 0 is the frequency of incident light. T is the ambient temperature, h is Planck's constant, c is the speed of light, k is Boltzmann's constant, and C is an arbitrary constant. The conversion process above was conducted via Multiwfn code [38]. The laser source wavelength used in our experimental platform is 532 nm, and room temperature was set as 298.15 K.

Comparison of Theoretical and Experimental Spectra
The experimentally measured and theoretically calculated spectra of DBDS, DBS, and BiBz are provided in Supplementary Material. It can be found from the Figures S4-S9 that the experimental Raman spectrum is in good agreement with the theoretical spectrum. Generally, the B3LYP hybrid functional combining with 6-311G** basis set overestimates the vibrational frequency in the high vibrational frequency zone but underestimates the frequency in the low frequency zone. This phenomenon is consistent with previous studies [43]. Dual scaling factor in the high frequency zone and low frequency zone is optional for agreement improvement.
In order to quantitatively describe the consistency between the theoretical spectrum and the experimental spectrum, a linear fitting operation is performed to each characteristic peak of DBDS. The abscissa of peak data is the experimental frequency, and the ordinate is the theoretical frequency. The peak data fitting line close to y = x indicates a better agreement between the experimental value and the theoretical value. As shown in Figure 6, the fitting lines of Raman spectrum and infrared spectrum both show good linearity. The goodness of fit reaches 0.99979 and 0.99981, respectively. A good agreement is revealed at B3LYP-D3(BJ)/6-311G** level as the slope of the fitted line is close to 1. The mode identification is shown in Tables 5 and 6. Most experimental Raman peaks can be identified to corresponding normal modes within a reasonable deviation. There are also some normal vibration modes which are not measured in the experimental spectrum due to the weak Raman activity or infrared activity.  Different conformations contribute to different experimental peaks. For example, in the both Raman and infrared experimental spectra, the peak at 2964 cm −1 is only generated by the #13 normal mode of DBDS-2, while the #13 normal mode of DBDS-1 contribute to the peak at 2908 cm −1 together with the #14 mode. The characteristic peaks generated by a single conformation are marked with parentheses in Tables 5 and 6.
In addition to the fundamental frequency band in the vibrational spectroscopy, there are many anharmonic phenomenon such as overtone, combination band, Fermi resonance, etc. Which contribute to several characteristic peaks in experiment spectrum that are not observed in the theoretical spectrum. For example, two peaks are observed in experimental Raman spectrum at 1495 cm −1 and 1526 cm −1 . While the adjacent peaks at 1463 cm −1 and 1584 cm −1 correspond to the #19, #20, #17, and #18 mode, respectively. There is no corresponding theoretical normal mode for these two peaks (1495 and 1526) assignment. They are actually a combination band of peaks at 1002 cm −1 and 481 cm −1 . The peak at 1002 cm −1 (#43 + #44) corresponds to in-plane bending vibration of benzene ring carbon skeleton, while the peak at 481 cm −1 (#67 + #68) is corresponding to the out-of-plane deformation vibration of the benzene ring. This kind of highly correlated vibration is prone to coupling and generating combined frequency band. Fermi resonance most often occurs between fundamental and overtone (combination) excitations, if they are nearly coincident in energy. The two peaks (1495 and 1526) are thus split under Fermi resonance with the action of strong fundamental frequency band (1600 cm −1 ).

Comparison of DBDS, DBS, and BiBz
DBDS, DBS, and Bibz have remarkable similarities in molecular structure. The only difference is the number of sulfur atoms between the two benzyl groups. The electronegativity of sulfur atoms is stronger than that of carbon atoms and hydrogen atoms. The sulfur atom will cause the redistribution of electron as discussed in Section 3.1, which changes the molecular vibration frequency. The spectra comparison of DBDS, DBS, and BiBz is shown in Figure 7. In order to reveal the similarity and difference of the vibrational behaviors of the three molecules, the six vibrational behaviors of the methylene groups is analyzed. The frequency distributions of these six vibrations are generally consistent. The order from largest to smallest is asymmetric stretching, symmetric stretching, scissoring, wagging, twisting, and rocking, according to the frequency magnitude.
The characteristic vibration frequency of six different vibration behaviors in the DBDS, DBS, and BiBz is shown in Table 7. The trend is plotted in Figure 8 and the frequency is averaged by two dual vibrations. It can be found that, for the stretching vibration located in the high frequency region, when the sulfur atom in the middle of the benzyl group is lost, the stretching vibration frequency tends to decrease, while the four bending vibration frequencies in the middle frequency region tend to increase.
The most interesting aspect of these graphs is the low rocking vibration value of BiBz-2, which shows a different trend from several other conformations. This is because the frequency difference between the two dual vibration is very large. BiBz-2 is the only molecule with C 2h symmetry among the five dominant conformations. The rocking vibration at 748 cm −1 and 964 cm −1 are synchronous rocking and asynchronous rocking, respectively. The synchronous rocking at 748 cm −1 keeps symmetry during the vibration based on the center plane of two benzene rings. The polariseability change rate of this vibration is zero, thus the Raman activity is zero. Exhibiting extremely low frequency.

Conclusions
In this paper, the vibration spectroscopy (Raman and infrared) and molecular vibrational behaviors of widely concerned molecules in sulfur corrosion phenomenon (Dibenzyl Disulfide, Dibenzyl Sulphide, and Bibenzyl) is detailedly analyzed. First, the proportion of different conformations are calculated according to Boltzmann distribution in relative Gibbs free energy via Molclus program. In order to obtain an accurate Gibbs free energy, the single-point energy is calculated in revDSD-PBEP86-D3(BJ)/may-cc-PVTZ and thermodynamic correction is calculated in the B3LYP-D3(BJ)/6-311G** level, respectively. Five dominant conformations of these three molecules were obtained, and their intramolecular interactions were detailedly discussed by NCI analysis.
PED based vibrational assignment of DBDS, DBS, and BiBz was performed via VEDA4 program, the normal vibration modes are assigned. Conformations weighted spectra is compared with experimentally measured spectra. It can be found that experimental spectra are in good agreement with the theoretical spectra in the B3LYP-D3(BJ)/6-311G** level with a frequency correction factor (0.9640). The characteristic peak in experimental spectra is also identified to corresponding modes within a reasonable deviation. The overtone due to the anharmonic effect is also pointed out. The divergence among these molecules is discussed. It can be found that, for the stretching vibration of methylene group located in the high frequency region, when the sulfur atom in the middle of the benzyl group lost, the stretching vibration frequency tends to decrease, while the four bending vibration frequencies in the middle frequency region tend to increase.
In summary, the results presented in this paper provide a theoretical and experimental guidance for understanding the vibrational spectra and molecular vibrational behaviors of these molecules. We also believe that this article will lay a solid foundation for the development of advanced corrosive resistant materials and realization of the corrosive state monitoring.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms23041958/s1. Author Contributions: Z.W.: conceptualization, methodology, software, validation, visualization, investigation, and writing-original draft preparation; R.S.: conceptualization, methodology, and writing-review and editing; J.W. and P.W.: conceptualization; Z.Z. and X.Z.: writing-review and editing; W.C.: supervision and project administration; and F.W.: funding acquisition. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the National Natural Science Foundation of China under grant number 51977017.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.