1 Introduction

Polyaromatic heterocycles (PAH) are among the pollutants found in soils in many industrialized countries, especially in sites that have been contaminated with coal–tar. Significant concentrations of thiophene, benzothiophene, and dibenzothiophene have been found in groundwater at coal–tar polluted zones [1]. The pollutant migration through the soil is controlled by diffusion and sorption on the surface of soil components, especially clay minerals. These minerals with small particle size and large surface area dominate the sorption of pollutants and organic matter in soils [2]. Besides, these polyaromatic hydrocarbons can be detected in atmospheric aerosols and phyllosilicates have been detected in dust-aerosols. Dust aerosols may be potential reactive surfaces for heterogeneous chemistry and a significant sink for many volatile organic compounds, as well as a vehicle for transportation of trace gases through the atmosphere. Although dust particles are very abundant in the atmospheres of Earth and Mars, comprising nearly 45% of global aerosol emissions in Earth [3], the role of mineral aerosols in atmospheric chemistry is currently not well understood. Therefore, it is interesting to study the adsorption process of thiophene derivatives on phyllosilicate surfaces at molecular level.

On the other hand, thiophenes are present in coal–tar and petroleum. Deep removal of organosulfur compounds from fuels is highly demanded, due to their poisoning effect on catalytic processes in the process of refining these fuels. Aromatic thiophene derivatives are especially difficult to remove. Different hydrodesulfuration processes used have the intrinsic problem that they reduce the octane index in fuels. An alternative is the adsorptive removal of organosulfur components. Many types of adsorbents, such as metal oxides, active carbon, clays, and zeolites [4], have been reported for the adsorptive removal of organosulfur compounds. The use of natural clay minerals for this application could represent an interesting environmental and cheap solution.

In addition, some thiophenes are important building blocks for materials that are of interest in nonlinear optics and that have applications in electronic devices [5]. The formation of nanostructured composite materials with these heterocycles and layered solids, like clays, can provide interesting physical, chemical, mechanical, and electronic properties [6].

A simple way of monitoring the presence of thiophenes in any media involves the use of infrared (IR) spectroscopy; hence, the determination of the spectroscopic properties of these compounds and their behavior upon adsorption processes is interesting. Previous authors have reported some variations in the IR spectra of thiophene [4, 7, 8] when they are adsorbed on silica, zeolite, and clay [9] surfaces.

In all these processes, the thiophene-phyllosilicate adsorption interactions are important. However, it is difficult to systematize data of these interactions experimentally, in part because the data come from different techniques and conditions. In this context, theoretical calculations can be a useful tool in understanding the experimental behavior. Quantum-mechanic calculations have been applied to the study of these molecules [10, 11]. Also atomistic calculations have been applied to the study of adsorption of organics on silicate surfaces, [12–16] and other minerals [17]. In this work, new theoretical approaches are used to study the molecular and electronic structures and spectroscopic properties of thiophene, benzothiophene, and dibenzothiophene (Fig. 1) and the effect of the adsorption of thiophene on a pyrophyllite surface on the spectroscopic properties.

Fig. 1
figure 1

Molecular structures of thiophene (a), benzothiophene (b), and dibenzothiophene (c)

Full size image

2 Models and computational methodology

2.1 Molecular calculations

Quantum-mechanic calculations of the thiophene derivatives were performed using the density functional theory (DFT) and Hartree–Fock methods. The electronic structure of these molecules was calculated using the hybrid functional BHandHLYP and the second-order Moeller–Plesset approximation at the all electron MP2 (full) level, both with a triple-ζ 6-311G** [18, 19] basis set with polarization functions on all atoms including H atoms. The Gaussian03 program package [20] was employed. All geometries were fully optimized at these levels using the Berny analytical gradient method. No geometry constraint was applied on the atoms of the molecules. Vibration frequencies of the normal modes were obtained from the analysis of the forces yielded from finite atomic displacements. All optimized structures showed only positive eigenvalues as corresponding to actual minima.

2.2 Periodic-model calculations

As a preliminary calculation for the adsorption of thiophene on mineral surface, we calculate the thiophene molecule as a periodical model. Then the thiophene was introduced in the center of a box with periodical boundary conditions with the lattice vectors: a = 10.5 Å, b = 9.1 Å, c = 20 Å, α = 90°, β = 100°, and γ = 90°. This lattice box is constant during the geometrical optimization (constant volume). This cell is large enough for avoiding intermolecular interactions with vicinal cells.

A periodic model of a phyllosilicate crystal structure was generated and used for our calculations. Pyrophyllite is a dioctahedral 2:1 phyllosilicate in which a sheet of octahedrally coordinated Al cations is sandwiched between two sheets of linked silica tetrahedra by bridging O atoms (named basal oxygens). The octahedral sheet is bonded to the tetrahedral ones by apical structural oxygen atoms, and the sixfold coordination of the aluminum atoms is completed by –OH groups. This triple-sheet forms a layer separated from another similar one by an interlayer space. The basal oxygens of the tetrahedral sheet form the surface of this layer. Experimental atomic coordinates and lattice cell parameters of pyrophyllite were taken for the crystal structure model [21]. All octahedral cations are Al and all tetrahedral cations are Si, and no interlayer cations were included. The H atom coordinates were included manually and optimized previously by means of interatomic empirical potential calculations [22]. This structure is in the trans-vacant crystal form, [23] and completely dry in our simulations in order to avoid too many variables for our systematic study. Before studying the adsorption process, we optimized completely the crystal lattice of the pyrophyllite. The structure was treated as triclinic (P1 symmetry) and all cell parameters of the lattice vectors, a, b, c, α, β, and γ were considered independent variables during this optimization.

For the study of the adsorption process, we use the surface corresponding to the 001 plane of the pyrophyllite, because this is the most stable surface and it corresponds to the largest specific surface that can be the external one or the interlayer one. In order to generate a surface along the 001 plane as a model of the external surface of this mineral, an additional vacuum space of 13 Å along the c axis was created in the crystal lattice previously optimized.

On the other hand, the maximum molecular size of the thiophene is 4.5 Å (from geometry optimized at BHandHLYP/6-311G**) and the minimum size in one unit cell of pyrophyllite is only 5.2 Å along the a lattice vector. Hence, for the purpose of the investigation of the interaction between the molecule and the mineral surface, a 2 × 1 × 1 supercell was generated for neglecting intermolecular interactions at the surface. This model shows a surface of 10.5 × 9.1 Å2, approximately 1 nm2 and the coverage density is one molecule per nm2 that represents the maximum coverage for a planar monolayer of thiophene on the mineral surface.

Energy calculations of the periodic crystal models were performed using the numerical atomic orbital (NAO) methodology implemented in the SIESTA program [24]. This is a method that scales linearly with the number of atoms in the simulation cell, and it is based on DFT and pseudopotentials. The generalized gradient approximation (GGA) and the Perdew–Burke–Ernzerhof (PBE) parameterization of the exchange-correlation functional were used. A uniform mesh with a plane-wave cut-off energy is used to represent the electron density, the local part of the pseudopotential, and the Hartree and exchange-correlation potentials. During the geometry optimizations a mesh-cut-off energy of 150 Ry was used and finally a total energy calculation with a mesh-cut-off energy of 500 Ry was performed according to previous work [25]. Core electrons have been replaced by norm-conserving pseudopotentials [26]. The pseudopotentials simulate the interaction between the valence electrons and the cores (nuclei plus core electrons) and neither core electrons nor core wave functions have to be included explicitly. With this approximation the valence wave functions are substituted by pseudo-wave functions that do not present strong oscillations in the core region. The basis sets are made of strictly localized NAOs. Their localization cut-off radii correspond to an energy shift of 270 meV. The basis set used in this work is double-Z polarized (DZP) following the perturbative polarization scheme. Calculations were restricted to two k-points of the Brillouin zone according to previous work [27]. In all structures, all atoms were relaxed by means of conjugated gradient minimizations at constant volume. The vibration frequencies of the mineral crystal lattice and surface were obtained with the VIBRA program included in the SIESTA package. Vibration frequencies of the normal modes were obtained from the analysis of the forces yielded from finite atomic displacements. The analysis of normal modes was performed with the Molekel program [28].

3 Results and discussion

3.1 Molecular structure of organics

In order to compare several quantum-mechanical methods, the thiophene molecule was optimized at the BHandHLYP/6-311G** and MP2/6-311G** levels with GAUSSIAN03, and with the PBE/DZP method using SIESTA. In all cases, a full molecular geometry optimization was achieved, where all molecular variables were named differently and allowed to vary independently. In Tables 1, 2 and 3, the main geometrical parameters of these optimized structures are compared with those of previous quantum-mechanical works [10, 11] and with experimental data [29].

Table 1 Main geometrical features (bond lengths (Å) and angles (°)) of calculated and experimental thiophene molecular structure
Full size table
Table 2 Main geometrical features (bond lengths (Å) and angles (°)) of calculated molecular structure of benzothiophene and dibenzothiophene (for all features see the Supplementary Material)
Full size table
Table 3 Mulliken net atomic charges (in electrons), dipole moment (in Debyes), and energy of HOMO and LUMO orbitals (in atomic units) of thiophene optimized at different calculation levels
Full size table

3.1.1 Molecular geometry

Experimental structure data are available for free thiophene [29] in the gas phase from microwave data, and for dibenzothiophene in the solid state from X-ray diffraction data [23]. Our thiophene calculations reproduce quite well the experimental C–S and C–H bond distances, and our results are closer to experimental data than those reported in previous calculations (Table 1). The calculations at the MP2/6-311G** level yield the geometry features that are closest to the experimental data. All bond distances calculated with SIESTA are longer than those based on 6-311G** level, while the main SIESTA angles are very close to experiment. The C–S, C2–C3, and C–H bond lengths are within the values of MP2/6-311G** and BHandHLYP/6-311G** levels. Similar tendencies and qualitative differences between calculation levels are observed in the benzothiophene molecular structure (Table 2; Supplementary Material). Taking into account the thiophene results, we can estimate that the experimental value of the C–S bond length in benzothiophene should be slightly longer than 1.730 Å, within the values of MP2/6-311G** and BHandHLYP/6-311G** levels, in the 1.730–1.737 Å range. Analogously, the experimental value of the C2–C3, and C–H bond lengths will be within the values of MP2/6-311G** and BHandHLYP/6-311G** calculations, the C2–C3 one being slightly shorter than 1.369 Å in the 1.369–1.342 Å range and the C2–H one slightly shorter than 1.083 Å in the 1.083–1.072 Å range. In thiophene and benzothiophene, the C2–H bond length at the α-position with respect to the S atom, is shorter than the rest of the C–H bond lengths, in agreement with experimental data [29]. The S–C9–C4 angle (111.41° at MP2/6-311G**) is smaller than the S–C2–C3 one (113.46°) due to the big value of the bicyclic C8–C9–S angle (127.13°). The other bicyclic angle is even bigger (C3–C4–C5 = 129.43°) forcing the low value of the C3–C4–C9 angle (111.63°). Similar differences between calculation levels are observed in the main geometry features of the dibenzothiophene molecule (Table 2 and Supplementary Material), where the calculations based on 6-311G** level yield geometrical features that are closer to experiment than the ones obtained in previous calculations [30]. Nevertheless, the calculations based on SIESTA and MP2/6-311G** levels yield the closest geometrical features to experimental data, especially for the C–S bond length in dibenzothiophene. In these compounds, the C–S bond length follows the following sequence: thiophene < benzothiophene < dibenzothiophene, due to the π electron delocalization on the aromatic ring.

Our calculations reproduce quite well the dibenzothiophene X-ray diffraction data [31], especially the C–S bond length. This result indicates that the intermolecular interactions in the dibenzothiophene crystal lattice are weak. Hence, the molecular structure of dibenzothiophene in the gas phase is similar to that in crystal form.

3.1.2 Electronic structure

The analysis of the charge distribution can be performed observing the relative values of the Mulliken net atomic charges in thiophene (Table 3), benzothiophene, and dibenzothiophene (Table 4 and Supplementary Material). In thiophene no significant differences are observed between results obtained with the BHandHLYP/6-311G** and MP2/6-311G** methods, although the MP2 calculations tend to yield slightly larger both negative and positive charges. Calculations using the 6-311G** basis set yield lower charges than previous calculations that used a smaller basis set (6-31G**), [10] indicating that bonds are less polarized. The S atom is positively charged and clearly represents the main nucleophilic attack centre. The C atoms in α positions with respect to the S atom have the largest negative charges and are, as expected, reception centers for an electrophilic attack [10]. In the same way, the H atoms of the α C atoms have a higher positive charge than those of the β C atoms.

Table 4 Main values of the Mulliken net atomic charges (in electrons), dipole moment (in Debyes), and energy of HOMO and LUMO orbitals (in atomic units) of benzothiophene and dibenzothiophene (for all values see the Supplementary Material)
Full size table

The total dipole moment calculated for thiophene is larger than the experimental value (0.55 D) [32] and considerably closer to experiment than in previous reported calculations [10]. The calculated (BHandHLYP/6-311G**) μ value is closer to experiment than the MP2/6-311G** one. With respect to the energy of the HOMO and LUMO orbitals, large differences occur between the results obtained with both methods: while the MP2 E HOMO value is more negative than BHandHLYP, the corresponding E LUMO value is much more positive, leading to an E HOMO–E LUMO difference significantly larger in MP2/6-311G** than in BHandHLYP/6-311G** (Table 3). This change is expected to affect the reactivity of this heterocycle with nucleophiles.

The introduction of a benzene ring in benzothiophene decreases slightly both the S atom positive charge and the heterocyclic C atoms negative charge, due to the possibility of π-electron delocalization. A similar effect is observed in dibenzothiophene. Nevertheless, the relative differences in the net atomic charges are similar in thiophene, benzothiophene, and dibenzothiophene; that is, addition of benzene rings does not change significantly the electronic structure of the heterocycle. This is in agreement with the results obtained by previous authors [10]. As in thiophene, the effect of the calculation level is not significant, in relative terms, for the description of the electronic structure in benzothiophene and dibenzothiophene. The dipole moment (μ) value calculated at MP2/6-311G** level is similar to the one reported in previous MP2/6-31 + G* calculations [10]. The μ value increases following the sequence: thiophene < benzothiophene < dibenzothiophene. In dibenzothiophene, the calculated μ value is larger than the experimental one [33], as found above in thiophene. In all cases the calculated dipole moment is larger with MP2/6-311G** than with BHandHLYP/6-311G**, and the difference increases also with the sequence: thiophene < benzothiophene < dibenzothiophene. Finally, the (E HOMO–E LUMO) gap decreases following the sequence: thiophene > benzothiophene > dibenzothiophene. This can be explained by the increase of the aromatic electron delocalization and hence the decrease of the energy difference between electronic levels increasing the absorption capacity of electromagnetic radiation.

3.1.3 Spectroscopic properties

The main vibration modes of the thiophene molecule, optimized at different levels, are reported in Table 5 together with the available experimental data [34]. In general, all calculated frequency values are slightly larger than the experimental values, the ones based on the SIESTA method being the closest to experimental data. For all modes, the frequencies calculated at the BHandHLYP/6-311G** level are larger than those found with MP2/6-311G**. Nevertheless, excellent linear relationships between calculated and experimental frequencies are obtained (Fig. 2) with scaling factors that are close to unity: 0.974 (correlation coefficient: R = 0.9998), 0.960 (R = 0.9998), and 0.938 (R = 0.9994) for the SIESTA, MP2/6-311G**, and BHandBHLYP/6-311G** calculations, respectively. In the particular 1,300–400 cm−1 range, the SIESTA calculations reproduce the experimental frequencies with a linear scaling factor of unity (exp. = 1.00 × calc., R = 0.9979).

Table 5 Calculated and experimental frequencies (cm−1) of the vibration normal modes of thiophene and thermodynamic properties (in cal/K mol) at 298.15 K and 1 atm
Full size table
Fig. 2
figure 2

Relationship between experimental and calculated (with Siesta) frequencies of thiophene

Full size image

The assignment of the thiophene vibration modes is consistent with the one proposed in experimental work, both in terms of frequencies and intensities [34]. In addition, our calculations allow us to complete previous assignments. In the high-frequency bands corresponding to stretching ν(CH) vibrations, our calculations can clearly distinguish between the different modes. The ν(CH) bands corresponding to C atoms in α positions with respect to S atom occur at higher frequencies than those of β C atoms. This is consistent with the shorter bond length and higher polarization of the Cα–H bond (Tables 1, 3). The symmetric and asymmetric C–H modes can also be distinguished in our calculations, though frequency differences are very small. The bands that appear in the 1,650–1,350 cm−1 range can be assigned to ν(CC), and they correspond to ring stretching, in agreement with the experimental νR bands. The 1,232 cm−1 band in SIESTA calculations corresponds to the experimentally assigned δR (ring) and involves bending atomic motions both in the ring and in C–H. We have tentatively assigned it to a coupling between δR (ring) and δ(CH). The bands in the 1,085–1,053 cm−1 range (in SIESTA calculations) should be assigned to δR modes. Although the experimental assignment for this group of bands was δ(CH), we observe that the participation of CC bonds in the atomic motions is also significant.

The most intense bands occur at 836 and 738 cm−1 (in SIESTA), and we have assigned them to the stretching ν(CS) and bending-out-plane γ(CH) modes, respectively, being consistent with the 839 and 712 cm−1 reported in the experimental data. The large intensity of the ν(CS) mode can be due to the high polarity of the C–S bond, as shown above in the net atomic charges (Table 3), and to the large variation in μ during this stretching mode. In the γ(CH) symmetric mode, all the C–H bonds change their orientation with respect to the aromatic ring simultaneously, thus producing a large change in the dipole moment that explains the very strong band intensity observed.

The benzothiophene and dibenzothiophene vibration mode frequencies are given in Tables 6 and 7, respectively. The ν(CH) bands of the heterocyclic CH bonds (C2–H, and C3–H) occur at higher frequencies than those of the benzene rings CH bonds (C5–H,···). Moreover, the heterocyclic CH bonds in the substituted thiophenes occur at frequencies that are almost 200 cm−1 larger than in thiophene. Nevertheless, the frequency differences between ν(C2–H) and ν(C3–H) are similar in thiophene and benzothiophene (40 cm−1). Symmetric modes always present slightly higher frequency than the asymmetric ones. Bands in the range 1,715–1,619 cm−1 in benzothiophene and 1,719–1,708 cm−1 in dibenzothiophene correspond to ν(CC) modes. They have larger frequencies in the benzene ring than in the thiophene ring (Table 5). In contrast with thiophene, the limit between ν(CC), δ(CH), and δR bands in benzothiophene and dibenzothiophene cannot be distinguished clearly. This fact is due to that the benzene δ(CH) modes occur at higher frequencies (1,550–1,518 cm−1 in benzothiophene) than in the heterocyclic ring (1,153 cm−1 in benzothiophene). The δR ring modes occur at frequencies up to 831 cm−1 in benzothiophene and 1,086 cm−1 in dibenzothiophene. At lower frequencies, the γR bands appear, but the frontier with δR and ν(CS) modes is not easy to distinguish clearly. In benzothiophene, bands in the 1,421–1,218 cm−1 range correspond to several coupled δ(CH) and νR modes. In dibenzothiophene, three frequency ranges can be distinguished: 1,684–1,672 cm−1, where the contribution of νR is bigger than that of δ(CH); the 1,555–1,515 cm−1 range where δ(CH) participates more than the νR mode; and the 1,389–1,333 cm−1 where both modes are coupled evenly.

Table 6 Frequencies of the vibration normal modes (cm−1) and thermodynamic properties (in cal/K mol) at 298.15 K and 1 atm of benzothiophene calculated at BHandHLYP/6-311G** level
Full size table
Table 7 Frequencies of the vibration normal modes (in cm−1) and thermodynamic properties (in cal/K mol at 298.15 K and 1 atm) of dibenzothiophene calculated at BHandHLYP/6-311G** level
Full size table

In Fig. 3, the calculated IR spectra of these molecules are described. In thiophene, we can distinguish clearly the ν(CS) and γR bands, whereas these cannot be distinguished in benzothiophene and dibenzothiophene. This fact can be explained due to the shorter bond length and higher polarity of the C–S bond in thiophene.

Fig. 3
figure 3

Calculated infrared spectra from BHandHLYP/6-311G** calculations of thiophene (a), benzothiophene (b), and dibenzothiophene (c)

Full size image

From the vibration data it is possible to calculate the thermodynamic properties of the three molecules considered. The entropy (S°) and heat capacity (C p°) calculated at 298.15 K and 1 atm reproduce well the experimental values for thiophene (Table 5) [34]. These properties have also been calculated in benzothiophene (Table 6) and dibenzothiophene (Table 7), though no experimental values are available for comparison. Nevertheless, the good agreement obtained in thiophene validates the theoretical values for the other molecules for future experimental data. Both properties increase in the following sequence: thiophene < benzothiophene < dibenzothiophene.

3.2 Adsorption of thiophene on phyllosilicate surface

3.2.1 Mineral surface and thiophene adsorption

The periodic model of pyrophyllite was fully optimized with the SIESTA program, and the following cell parameters were obtained: a = 5.25 Å, b = 9.10 Å, c = 9.32 Å, α = 90.1°, β = 100.8°, γ = 89.9°. These values are close to the experimental values (a = 5.16 Å, b = 8.97 Å, c = 9.35 Å, α = 91.2°, β = 100.5°, γ = 89.6°) [21]. The average value of the Si–O bond length is 1.66 Å, which is slightly longer than the experimental value of 1.62 Å. The average value of the Al–O and Al–OH bond lengths are 1.95 and 1.92 Å, respectively, and reproduce well the experimental values [21]. The OH bond length is 0.979 Å with an inclination angle ρ of 30–32° that is consistent with experimental data. From this optimized structure, a 2 × 1 × 1 supercell was generated and a surface on the 001 plane was created with a c parameter value of 20 Å.

Next, the thiophene molecule previously optimized at BHandHLYP/6-311G** level was added at a distance of 3 Å to the pyrophyllite surface (layer of basal oxygens) in a parallel position with respect to the (001) surface, as suggested in previous experimental work [4]. The whole complex, surface + adsorbate, was then fully optimized at constant volume with SIESTA (lattice cell parameters: a = 10.5 Å, b = 9.10 Å, c = 20 Å, α = 90.1°, β = 100.8°, γ = 89.9°). This allows to the thiophene to search the minimum energy configuration with respect to the surface optimizing at the same time the molecular structure in the new situation. We found that, in the optimized adsorption complex, the thiophene molecule adopts an orientation parallel to the mineral surface at an approximate distance of 2.7–3.3 Å, and situated rather in the centre of one tetrahedral ring (in Fig. 4b the tetrahedral ring is highlighted for showing the relative position of the thiophene molecule). In general the S and H atoms, that are positively charged, lie above the basal O atoms of the pyrophyllite surface (Fig. 4). The S atom occupies a position above a Si atom and at 3.164, 3.390, and 3.660 Å to the three O atoms joined to this Si atom. The H atoms are in a balanced position where each one is at a non-bonding distance smaller than 4 Å to at least three O atoms. No significant variation in the molecular structure of thiophene is observed upon adsorption, except for the S–C bond length, which is slightly longer than in the gas phase (1.725 Å, Table 1), and a certain asymmetry is observed in the adsorbate molecule: d(S–C2) = 1.731 Å, d(S–C5) = 1.728 Å. The latter is due to the different local environment in the adsorption complex, where the shortest non-bonding distances to the surface are: d(C2H···O) = 3.277 Å, d(C3H···O) = 2.838 Å, d(C4H···O) = 2.925 Å, d(C5H···O) = 2.720 Å. The d(C–H) bond distances corresponding to the C atoms in the α position with respect to S are slightly longer (1.096 Å) than in the gas phase (1.094 Å), due to the interaction of these H atoms with the mineral surface O atoms.

Fig. 4
figure 4

Adsorption complex (pyrophyllite + thiophene) optimized with SIESTA, views from the (100) (a) and (001) (b) planes of the crystal lattice of the mineral (the orientation of the coordinate system is showed in the low-right part of each structure). The H atoms are highlighted as white balls, except the one below the adsorbate that is colored in blue. In (a) the O atoms that are close to the blue H atom are highlighted as red balls, the basal (small red ball) and the apical (big red ball) O atoms. The Al, Si, C, S, and O atoms are in pink, orange, gray, yellow, and red colors, respectively

Full size image

The interactions in the adsorption of this aromatic compound on the hydrophobic surface of pyrophyllite are weak and the adsorption distances are large as they correspond to mainly electrostatic and van der Waals interactions. Hence, the potential energy surface of this adsorption is more or less flat. In order to explore other possible minimum energy configurations, different orientations of thiophene molecule to the one described above (M1, Fig. 4) were probed (Fig. 5): M2, after starting with the thiophene at 5 Å parallel to the surface, the optimization yielded M2 where the thiophene is planar with a slight rotation of the aromatic plane with respect to the pyrophyllite surface and the aromatic ring is slightly displaced from the tetrahedral vacancy with respect to M1; M3 where the thiophene ring is also planar with respect to the surface but the ring is over the bridging O atoms between two tetrahedral cavities; M4 where the thiophene ring is in a perpendicular disposition with respect to the surface and the H atoms of the C3 and C4 atoms are oriented towards the basal O atoms of the surface; and M5 like M4 but the S atom is oriented to the surface. All configurations were fully optimized at constant volume and the configuration of minimum energy was M1, the rest being 1.7 (M2), 2.2 (M3), 7.0 (M4), and 7.9 (M5) kcal/mol less stable than M1.

Fig. 5
figure 5

Alternative adsorption complexes (pyrophyllite + thiophene) explored and optimized with SIESTA view from the (100) plane of the crystal lattice of the mineral. The Al, Si, H (highlighted in balls), C (in balls), S (in balls), and O atoms are in pink, orange, gray, white, yellow, and red colors, respectively

Full size image

The limitations of the DFT-PBE approximation for describing weak interactions, as in our case, are already known and several new functionals are currently appearing to overcome this problem [35, 36]. However, in our system the adsorption interactions are stronger than in pure aromatic hydrocarbons, like benzene or graphene. The presence of the heteroatom in thiophene produces an electron polarization and electrostatic interactions appear in the adsorption on the silicate surface. The adsorption energy for the most stable adsorption complex M1 is −15.75 kcal/mol that is slightly higher, though in the same order, than previous DFT calculations on adsorption of dioxins on pyrophyllite (−11.5 kcal/mol) [16]. Nevertheless, our calculated geometry of the adsorbate on the mineral surface is consistent with experimental data [4] and previous calculations on adsorption of other aromatics, like benzene [15] or dioxin [16] on phyllosilicates. A further discussion on the adsorption energy and calculation levels is beyond the scope of this work, especially for large systems like ours. Nevertheless, our approach yields quite better results for spectroscopic study, closer to experimental frequency values, in our systems than other calculations based on Hartree–Fock [37] or empirical potentials [22] in phyllosilicates.

3.2.2 Spectroscopic properties

In the thiophene–pyrophyllite adsorption complex, the ν(OH) modes of pyrophyllite can be easily identified. Their frequencies fall in the 3,676–3,666 cm−1 range, in good agreement with the experimental value [38] for free pyrophyllite (ν(OH) = 3,675 cm−1). However, the ν(OH) mode corresponding to the OH group that is just below the adsorbed thiophene molecule occurs at a slightly lower frequency (3,630 cm−1) (in Fig. 4, the H atom of this OH group is in blue colour). This can be explained by observing the geometry of this particular OH group. The presence of the thiophene molecule changes slightly the OH bond length (0.979 Å instead of 0.978 Å of the other OH groups), as well as the inclination angle of the OH bond with respect to the 001 plane, ρ(OH) (Fig. 4a), (32.8° instead of 30.7°). This fact can be due to the weak electrostatic interactions between the positively charged H atom and the π-electron ring of the thiophene that is just above the H atom (Fig. 4b). These geometric changes produce a stronger hydrogen bridge with the tetrahedral basal O atom, which is in front of it at 2.26 Å (in Fig. 4a it is as a small red ball). Instead, the corresponding non-bonding H···O distance for the other OH groups is 2.34–2.39 Å. The new orientation causes a slight increase in the non-bonding H···O distances with the nearest apical oxygens (in Fig. 4a they are represented as medium-size red balls) (2.90–2.91 Å instead of 2.87–2.88 Å), and the interaction produces a shift towards lower frequencies on the ν(OH) modes with respect to the rest of OH groups. This phenomenon is observed also in mica and illites, where the presence of IVAl3+ in the tetrahedral sheet increases the strength of the hydrogen bridge and consequently decreases also the ν(OH) frequency [39, 40].

The thiophene ν(CH) modes of the molecule upon adsorption on pyrophyllite occur at slightly lower frequencies (Table 8) than in the gas phase (Table 5). This is due to the weak adsorption interactions. This was also observed experimentally in thiophene adsorbed on SiO2 [4], where these modes appear also at slightly lower frequencies than in the gas phase [34]. In the adsorbed molecule it is difficult to distinguish between the ν(CH) modes of each C atom because they are coupled to each other due to the adsorption interactions. The ν(C2-H) mode is the only one that can be identified separately at a slightly higher frequency, due to the fact that this C atom is a little further from the surface (see above) than the other C atoms, and thus the effect of adsorption is smaller. In experimental spectroscopic studies of thiophene adsorption, this shift to lower frequencies effect in the ν(C2H) mode can be observed best in montmorillonite[9] and even better in zeolites[4], because the adsorption interactions are stronger than in our case. A similar shift effect is observed when the thiophene molecule is in the liquid and crystalline state with respect to gas phase[41] because the intermolecular interactions are stronger in the crystal form. In the ν(CH) modes of the rest of the C atoms, this effect cannot be detected clearly because it is within the frequency range zone.

Table 8 Calculated (with SIESTA) and experimental frequencies of the vibration normal modes of the adsorption complex of thiophene adsorbed on pyrophyllite surface (cm−1)
Full size table

On the other hand, the ν(CC) modes of the thiophene ring occur at frequencies similar to those in the gas phase. A similar phenomenon was observed experimentally in thiophene adsorbed on SiO2, although these bands are spread on a wider range. In crystalline and liquid thiophene, a slight shift to lower frequencies (with respect to gas phase) was observed due also to the weakness of intermolecular interactions. In thiophene adsorbed on zeolites this interaction is stronger and a lower-frequency shift in this mode was observed experimentally upon adsorption (from 1,408 to 1,396 cm−1) [7].

The δ(CH) modes occur at slightly lower frequencies (1,074–1,067 cm−1) than in the gas phase (1,085–1,084 cm−1), as also observed in thiophene in the solid state [41]. Nevertheless, the ν(SiO) modes are coupled with δ(CH) and this conclusion can be considered only as tentative. The δ(CC), ν(CS), and γ(CH) modes are coupled with lattice vibrations.

The pyrophyllite bending-in-plane δ(OH) modes are coupled with lattice vibrations and appear at higher frequency than the experimental value reported in smectite/illites, due to the presence, in these minerals, of cation substitutions that do not occur in pyrophyllite [42]. In the pyrophyllite crystal lattice, the ν(Si–O) modes can be distinguished in the 1,099–983 cm−1 range, though sometimes they are coupled with thiophene ring modes. These values are consistent with experimental values in silicates in general, ν(Si–O) = 1,070–970 cm−1 [43]. Our calculations can also distinguish the ν(Si–O) modes from the apical O atoms, ν(SiOAl) vibrations.

4 Conclusions

This work has explored several theoretical methods for describing the molecular and electronic structure, spectroscopic, and thermodynamic properties of PAH. Our calculation levels are higher than those reported in previous work, and we are able to reproduce the available experimental values. In the Hartree–Fock approximation, a triple ζ basis set and polarization functions augmented with MP2 method and hybrid HF/DFT methods like BHandHLYP are necessary to obtain molecular structure and electronic properties close to experimental data. The MP2/6-311G** calculations yield a slightly better agreement with experimental geometrical features than BHandHLYP/6-311G** ones, whereas electronic properties (µ) are better described by BHandHLYP/6-311G**. Experimental thermodynamic properties, such as, entropy (S°) and heat capacity (Cp°) are well reproduced with both methods. The agreement between our calculated results and experiment, in terms of geometric, electronic, and thermodynamic properties, validates the application of these theoretical methods employed for studying and predicting these properties in PAH systems for which no experimental data is available.

The close similarity between the molecular structure of dibenzothiophene in the gas phase and the one in the crystalline state is remarkable. It indicates that the intermolecular interactions in the crystal lattice are weak and do not alter the molecular structure. This can be an interesting starting point to explore the molecular packing in the precipitation of asphaltenes, where PAH are present, which is an important problem in oil-extractions.

Using the DFT/GGA approximation, the PBE functional with pseudopotentials and NAOs, as implemented in the SIESTA program, reproduces quite well the experimental data, especially the spectroscopic properties. Although these calculations do not improve the geometric features obtained in the above BHandHLYP/6-311G** and MP2/6-311G** calculations, the DFT/GGA/PBE calculations reproduce very well the IR frequencies of the main normal modes of PAH. The excellent experimental spectroscopic data that are available for the thiophene molecule has allowed us to compare of normal vibrational modes with our calculations, and a good agreement is found. We find linear relationships between experimental and theoretical frequency values with scale factors in the range of 0.94 – 1.00, depending on the calculation level. Our theoretical results should be useful in spectroscopic studies for the accurate assignment of IR-bands, and they can be extended to other PAH for which the experimental assignments are more difficult.

A good advantage of the SIESTA calculations is that periodical systems can be studied. This point is important in PAH adsorption on solids studies, where the large molecular size of PAHs forces one to use large surfaces of supercells of crystalline solids. The calculation method is especially useful for phyllosilicates since it reproduces quite well the crystal lattice and surface of the mineral. However, there is a hot debate about the limitations of DFT and HF approximations for describing properly the weak Van der Waals interactions that exist in many adsorption processes, especially for absolute adsorption energy values, and several new methodologies are appearing in order to describe these weak interactions. It is not feasible for us to evaluate the inclusion of dispersion interactions on such large systems with current computational capacity. However, further calculations related with this point are the focus of current work.

Nevertheless, our good results in geometric and spectroscopic studies, make SIESTA a useful method for spectroscopic studies in adsorption of PAH’s on phyllosilicate surface. Our calculations show that adsorption on pyrophyllite surface produces a low-frequency shift of the ν(CH), δ(CH) bands and, in a smaller scale, in the thiophene ν(CC) modes, with respect to the gas phase. This effect is smaller than the one observed in the case of zeolites where the adsorption interactions are stronger because of the existence of cations and high electric gradients. Taking into account experimental studies of adsorption of thiophenes on more active surfaces, there is a direct relationship between the strength of the adsorption interactions and the low-frequency shift of the ν(CH), δ(CH) and ν(CC) modes with respect to the gas phase. No significant effect on the frequencies of pyrophyllite was observed, except for the OH groups that are closest to the adsorption sites, where a low-frequency shift of ν(OH) is observed. This shift is due to the adsorption interaction of the adsorbate, interactions of H atom of the OH group with the surrounding tetrahedral O atoms and the interaction of this H atom with the π-electron ring of thiophene. This result allows us to explain at atomic scale the changes observed in the infrared spectra and it can be very useful for further experimental work in the detection of adsorption of thiophenes on silicate surfaces by spectroscopic techniques.

This work is an initial step in the exploration of the adsorption of aromatics to the phyllosilicate surfaces. Further work with different conditions and the presence of other molecules (water, cations, etc..) will be performed.