Discerning Inter- and Intramolecular Vibrations of Sulfur Polyaromatic

Sep 2, 2017 - Institut des Sciences Analytiques et de Physico-Chimie pour l'environnement et les Matériaux, Équipe Chimie Physique, UMR 5254, Universi...
0 downloads 0 Views 5MB Size
Subscriber access provided by University of Glasgow Library

Article

Discerning Inter- and Intramolecular Vibrations of Sulfur Polyaromatic Compounds Patricia Guevara Level, Hugo Santos-Silva, Faustine Spillebout, Kirk H. Michaelian, John M. Shaw, Isabelle Baraille, and Didier Bégué J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b03903 • Publication Date (Web): 02 Sep 2017 Downloaded from http://pubs.acs.org on September 3, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Discerning Inter- And Intramolecular Vibrations Of Sulfur Polyaromatic Compounds Guevara Level, P.,† Santos Silva, H.,† Spillebout, F.,‡ Michaelian, K.H.,¶ Shaw, J.M.,‡ Baraille, I.,† and B´egu´e, D.⇤,† †Institut des Sciences Analytiques et de Physico-Chimie pour l’environnement et les ´ Mat´eriaux, Equipe Chimie Physique, UMR 5254, Universit´e de Pau et des Pays de l’Adour, Helioparc, 2 Avenue du President Angot, 64053 Pau Cedex 9, France ‡Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta, Canada, T6G 2G6 ¶Natural Resorces Canada, CanmetENERGY, One Oil Patch Drive, Devon, Alberta T9G 1A8, Canada E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Thiophenes are an important class of molecules in fields as diverse as petrochemistry, molecular electronics, and optoelectronics. Thiophenic submolecular motifs are thought to play a role in molecular association and nano-aggregation phenomena in both pure materials and natural and synthetic mixtures. Vibrational (infrared and Raman) spectroscopy provides the means to characterize these species. In this work far-infrared photoacoustic and low-frequency Raman spectra of a series of polycyclic aromatic hydrocarbons containing sulfur have been measured and interpreted using DFT calculations based on a perturbational-variational method coupled with potential truncation. The approach and outcomes illustrate how inter- and intramolecular vibrations for thiophenic systems in single and multicomponent mixtures can be discriminated. This work o↵ers the perspective to search the inter- and intra-molecular signatures of the main submolecular motifs and heteroelements postulated as being present in the asphaltenes.

Introduction Thiophenes are an important class of molecules in various domains of chemistry including petrochemistry, 1–3 molecular electronics, 4–7 and optoelectronics. 8 In particular, they are of great interest in all that concerns oil upgrading. 9–11 Their presence as aromatic fragments in asphaltenes has been revealed by methods such as thermolysis analysis, 12 mass spectroscopy, 13,14 and NMR and IR spectroscopies. 15,16 Importantly, thiophenes are hypothesized to play a key role in the aggregation of asphaltenes 17,18 which is of uppermost importance to the oil extraction process. Asphaltenes are known to cause the blocking of pipelines and refinery equipment.

The molecular structure of asphaltenes is difficult to determine since they are formed by thousands of molecules that tend to stick together in solution. Over the past six decades, the 2

ACS Paragon Plus Environment

Page 2 of 48

Page 3 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

composition, molecular structure and colloidal state of asphaltenes in crude oil and bitumen have been considered a central enigma of petroleum chemistry. 19 The noncovalent interactions responsible for molecular recognition among components determine the microstructure of the aggregates. Two contrasting molecular architectures have been proposed to describe the dominant molecular motifs in asphaltenes: the island model in which asphaltene molecules contain one polyaromatic core with pendant aliphatic chains; and the archipelago model in which several aromatic motifs are bridged together via aliphatic chains. 20,21 The interactions found in asphaltene aggregates are ruled by van der Waals forces, hydrogen bonds and electrostatic interactions among heteroatoms. This is the case regardless of the model to which the molecules belong (archipelago or island models). Porphyrins also play major roles. 22 The identification of interactions in asphaltenes, the understanding of how they influence molecular arrangements, and finally the experimental method by which they are identified, are major issues in production, transport, and refinery process design and operation. 23 For example, understanding the impact of inter- and intramolecular interactions occurring between thiophene-like molecules is a missing and necessary refinement to asphaltene aggregation models, 20,24,25 that will lead to designing better processes for cracking asphaltene aggregates and facilitate refining of thiophene-rich crude oils.

Two types of interactions coexist among aromatic molecules containing sulfur but ⇡-⇡stacking between aromatic rings, and ⇡-S between an aromatic ring and a sulfur atom are the principal interactions. The interaction between aromatic rings and sulfur atoms has been shown by Morgan et al 26 to play an important role in stabilizing the folded conformations of 8 di↵erent proteins. Tauer et al 27 studied the H2 S-benzene dimer system in order to visualize molecular structures resulting from ⇡-S interactions, and to establish an appropriate computational method capable of describing these interactions. The interactions driving the spatial arrangement of thiophene molecules have been studied theoretically by Tsuzuki et al 28 who showed that the major interaction among thiophene molecules was dispersive,

3

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

whereas the electrostatic interactions were responsible for the orientation of the molecules.

In our previous work, 29 we showed that interactions responsible for ⇡-⇡ stacking for acene molecules could be visualized in infrared and Raman spectra. The frequencies related to these interactions were found between 15 and 100 cm 1 . This range of frequencies is difficult to measure experimentally using conventional methods but with the development of terahertz spectroscopy, including synchrotron-based photacoustic (PA) IR measurements, it is now possible to identify absorption bands specific to intermolecular interactions that occur at very low frequencies and to evaluate their importance.

In the current study we (1) identify the presence of thiophene-type molecules in complex systems, (2) identify ⇡-stacking interactions when they are predominant in a complex system by means of vibrational spectroscopy experiments and high-precision quantum mechanical calculations, and (3) develop a new variation-perturbation strategy adapted to the resolution of the vibrational Schr¨odinger equation at very low wavenumbers for large-size systems (bigger than 15 atoms) and their dimers. These results are used to characterize and identify IR and Raman signatures of vibration modes found at very low wavenumbers that are characteristic of association, e.g. dimerization. More specifically, infrared and Raman spectra of benzothiophene 1a, dibenzothiophene 1b, 4.6-dimethyldibenzothiophene 1c and 4-methyldibenzothiophene 1d are reported (Figure 1). The observed bands are assigned using high level calculations incorporating anharmonic electrical and mechanical approximations. Fundamental bands previously assigned by Bree et al 30 are confirmed, and additional combination bands are identified. The approach and outcomes illustrate how inter- and intramolecular vibrations in single- and multi-component mixtures can be discriminated.

4

ACS Paragon Plus Environment

Page 4 of 48

Page 5 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 1: Molecular structures studied: from left to right: benzothiophene 1a, dibenzothiophene 1b, 4,6-dimethyldibenzothiophene 1c and 4-methyldibenzothiophene 1d

Experimental Spectra Dibenzothiophene 1b, 4.6-dimethyldibenzothiophene 1c and 4-methyldibenzothiophene 1d were obtained from commercial sources at purity levels of at least 98 per cent and analyzed as received. In some cases, samples were ground manually prior to acquisition of their spectra. Raman spectra

Raman spectra were obtained using two spectrometers at CanmetENERGY (Devon, Alberta). In one set of experiments, a Renishaw InVia microscope-based system was coupled to a Spectra Physics 125 He-Ne laser, which provided excitation at 633 nm. Several milligrams of each solid were examined with the use of a 50⇥microscope objective. Spectra were acquired with a spherically focused incident beam; typical resolution was better than 2 cm 1 . Total accumulation time was 5 min for each spectrum. The second group of experiments was performed with a Bruker IFS 88 FT spectrometer and FRA 106 Raman accessory. Excitation of the FT-Raman spectra was e↵ected with a 1064-nm Nd:YAG laser at power levels between 90 and 260 mW. Ten 32-scan spectra with a resolution of 8 cm

1

were averaged for

each sample.

Far-infrared spectra PA far-infrared spectra were acquired at a resolution of 6 cm

5

ACS Paragon Plus Environment

1

using two di↵erent Bruker

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 48

IFS 66v/S Fourier transform infrared (FT-IR) spectrometers. The first instrument, at CanmetENERGY, was operated in both rapid- and step-scan modes. The He-Ne laser modulation frequency was 1.6 kHz for the rapid-scan measurements. Large amplitude phase modulation 31 was employed in the step-scan experiments, with a Signal Recovery lock-in amplifier being utilized for demodulation. The second spectrometer was located at the Canadian Light Source (Saskatoon, Saskatchewan). Rapid-scan spectra were acquired under conditions similar to those just mentioned. A thermal (globar) source was employed with this instrument.

Standard MTEC 300 gas-microphone PA accessories were used with both spectrometers; the sample cells were fitted with polyethylene windows. Helium was employed as the carrier gas. Multilayer mylar beamsplitters were installed in the interferometers. Spectra were recorded between about 50 and 700 cm 1 , due to limits imposed by the optical materials. Carbon black powder or an MTEC carbon black reference yielded reference spectra that were used to correct sample spectra for the wavenumber-dependent response of each instrument.

Mid-infrared spectra PA mid-infrared spectra were recorded at a resolution of 6 cm

1

using the Bruker IFS 88

spectrometer mentioned above. The laser frequency was either 1.6 or 2.2 kHz. An MTEC 200 PA cell was utilized. Nitrogen was used to purge the spectrometer and as a carrier gas. A KBr window sealed the PA cell, and a Ge/KBr beamsplitter was used in the spectrometer. Ten 32-scan spectra were averaged for each sample. Spectra were acquired from about 400 to 4000 cm 1 . Carbon black powder was used to obtain reference spectra in these experiments.

6

ACS Paragon Plus Environment

Page 7 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Calculation method for thiophene-type molecules and their dimers Structural details The molecular structures of thiophene-like molecules and their respective dimers were calculated using Density Functional Theory (DFT) with the !B97X-D exchange-correlation functional and the 6-311++G** basis set, implemented in the Gaussian 09 suite of programs. 32 The !B97X-D functional developed by Chai et al 33 is based on optimized longrange corrected hybrid density functionals, which employ 100% Hartree-Fock (HF) exchange for long-range electron-electron interaction correction. This method, tested by Salzner et al 34 for ⇡-conjugated oligomers, was also successfully benchmarked for interacting polycondensed aromatic molecules, 29 thus proving it overcomes the failure of standard DFT calculation approaches to describe long-distance interactions. The 6-311++G** basis set is of triple-zeta quality for the valence electrons. Among the basis sets tested by Wiberg, 35 including the Dunning correlation basis sets such as aug-cc-pVTZ and aug-cc-pVQZ, the 6-311++G** basis set gave accurate geometries and frequencies at a relatively small computational cost. The calculations use planar ‘C2 v symmetry’ for dibenzothiophene, 4,6-dimethyldibenzothiophene, and Cs for benzothiophene and 4-methyldibenzothiophene. The geometrical parameters of the molecules are benchmarked with experimental data when available.

For dimers, as an alternative to the supermolecular treatment, an approach based on symmetry-adapted perturbation theory (SAPT) 36 which utilizes the description of the interacting monomers in terms of Kohn-Sham (KS) orbitals was used in this work. The DFT-based SAPT approach, 37 denoted as SAPT-DFT in this article, is exact for all major components of the interaction energy (asymptotically for exchange interactions) in the sense that these energies would be exact if the DFT description of the monomers were exact. It provides very accurate results for small dimers such as (C6 H6 )2 where results for 7

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 48

the much more computationally expensive CCSD(T) method are well approximated. 38 The SAPT-DFT calculations were performed using the MOLPRO2012 package. 39 The PBE0 functional 40 was used as the DFT functional in SAPT, using the aug-VTZ basis set. All computational outcomes are reported in the SI.

Resolution of the vibrational Schr¨ odinger equation (RVSE) Position of the problem Without any consideration of molecular size, the theoretical process generally developed to solve the Vibrational Schr¨odinger Equation in the mechanical anharmonic hypothesis requires two restricted steps: the construction of the potential energy surface (PES), and the resolution of the vibrational equation (RVE) in order to obtain the energy levels and consequently the fundamental, overtone and harmonic wavenumbers describing the vibrational spectra. Calculations of the intensities in the electrical anharmonic hypothesis generally complete the analysis of the spectra.

Owing to the size of the systems studied in this work, we wondered about the precision we could get on the PES that include hundreds of thousands of terms. The search for an optimal analytical expression of the PES is an ongoing topic of discussion, even in cases of low-dimensional systems (3640 terms for 6 atoms; 92752 terms for 12 atoms ... without any consideration of symmetry). The skills and experience of some research groups 41 42 show that the PES of large-size systems can be developed by expressing the potential function in a simple way, i.e., a Taylor series expansion using curvilinear displacement coordinates. These series were generally truncated at fourth order. Quadratic, cubic, and quartic force constants were obtained by fitting the electronic energy data generally calculated by DFT methods for various nuclear configurations close to the respective optimized geometry.

8

ACS Paragon Plus Environment

Page 9 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Of course, vibrational spectroscopy modeling outcomes are directly dependent on the quality of the potential energy surface (PES) but also depend on the mathematical method used to solve the Schr¨odinger vibrational equation. It is generally assumed that variational or VCC (Vibrational Coupled Cluster) approaches take into account the main couplings between the 3N

6 vibrational modes of a given system. However, the use of such methods

is restricted to low-dimension systems (including less than 10 atoms) or highly symmetric systems. Among the accepted variational approaches, those of the conventional VCI type (SDT. . .) should first be cited. 43–45 In these approaches, all the excited vibrational states (single, double, triple. . . excitations) of a target active space (AS) are generated to form a diagonalisation space (DS) obtained with no particular selection criterion (except for the maximum degree of excitation). This approach is, to some extent, the most accurate since it ensures that no state is omitted in the variational process that contributes to description of the desired AS. On the other hand, it also generates DSs of huge dimensions that may be technically unusable (partly because the matrix to be diagonalised must be stored in memory) once the studied system contains more than ten atoms (these methods have been nowadays presented with more than 1010 Hartree products, see refs 46–48 for details on recent CI strategies evolutions and performences). Among the proposed solutions to avoid the diagonalisation of very large matrices, one should cite the successes obtained by the iterative methods such as the Lanczos algorithm 49–52 to calculate eigenvalues of a large, sparse matrix. However, these methods do not permit operation on larger systems at this time. Despite the ubiquity of iterative methods reported in the literature in the last decade, developments related to explicit processing of the conventional CI matrix have not been abandoned and continue to be active. Carrington has recently published the most accurate example that illustrates these recent developments for the acetonitrile (CH3 CN) molecule. 53 In this work, DSs of several hundred thousand terms are used. Constrained by these limitations, this researcher has chosen to add a selective criterion of the DS to his variational approach by applying a pruning approach. 53 Despite this restriction, the size of the DS remains very

9

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 48

large. Recently we have shown that this space can be reduced even more, but our work also highlights the fact that this reduction does not allow the treatment of very much larger systems if one requires exact convergence of the variational problem. 54 In conclusion, even if one combines the pruning approaches to limit the expansion of the conventional VCI approaches (SDT. . . ) and the adaptative methods of the A-VCI, 55 or VCC, VCI, . . . , types, one can not ensure, currently, the possibility of studying systems containing more than 15 atoms by variational methods. Only the use of less restrictive variational methods like the VSCF type allow the attainment of such dimensions. 56 One should, nevertheless, note that new routes like those introduced by Carrington on tensors 57 or those developed in our recent work 55 are likely to tackle this challenge in the future. As these approaches have not yet been finalised, other ‘less precise’ solutions should be considered for dealing with medium-sized systems. Conversely, vibrational studies of larger systems are now common since the perturbative VPT2 approach is implemented in most quantum chemistry codes (Gaussian for example 58 ). To improve the description of perturbations, other combined approaches that attempt to take advantage of both variational and perturbational approaches are used. These are known as variation-perturbation methods. Within these approaches, only a few vibrational states in DS strongly coupling the AS state(s) for which the vibrational characteristics are to be determined are selected by a procedure of interactively enriching the diagonalization space while intentionally controlling (see ref 47 for a brief recall onto the VPT approaches ) the dimension of the space to be diagonalized. Finally, this approach amounts to the development of a VPT2 method in which it is possible to explicitly obtain a very accurate description of the main resonances, interactions, etc. by diagonalising DSs with only a few tens of thousands of terms.

Background, notations and definitions in relation to this method are reported in the next paragraph.

10

ACS Paragon Plus Environment

Page 11 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Theoretical Approach In vibrational structure theory, the standard choice of harmonic oscillator (HO) functions is widely used. We made this choice because, for the Taylor expansion approach, the one-mode operators are simply normal coordinates in powers up to the degree of the polynomial expansion. 55 The integrals required for the implementation of CI using a Hamiltonian written in the sum-of-product (SOP) form are one-mode integrals. In any case, there are only a few integrals per mode; the integral evaluation is one-dimensional and can thus be done simply and fast. Nevertheless, it is important to note that while HOs are an obvious choice for the expansion of the one-mode orbitals of harmonic vibrational modes, this is not true for general anharmonic molecular motions. In such cases the use of HO basis sets can cause instabilities in the variational optimization of the vibrational wave function, due to collapses in the configurational regions not properly described by the semi-global representation of the molecular PES. No variational collapses were detected by our code for the present study of the lowest frequency modes in the spectra, but this is a risk that must be kept in mind when working with HO-based sets. Another point that should be kept in mind when combining both HO and quartic PES conditions concerns the average wavenumber precision that can be expected. It is generally accepted (apart from any variational collapse) that the precision obtained is on the order of 0.5 to 1% with respect to the experimental data. 55,59–61

Both the definitions of all these terms and the principe of the method used in this work are now summarized below.

Let us consider the dimensionless normal coordinates q = (q1 , q2 , q3 , . . . qN ), where N is the number of vibrational degrees of freedom of the system. In the Born-Oppenheimer approximation, the vibrational Hamiltonian H(q) is the sum of the harmonic operator H0

11

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 48

and the anharmonic part V and is written as H(q) = H0 + V(q) ✓ ◆ N S N X X Y !i @2 2 = + qi + Ks qisi , 2 2 @qi i=1 i=1

(1)

ksk1 =3

where !i is the harmonic frequency (in cm 1 ) associated with the qi coordinate, ksk1 is the sum of elements of multi-index s = (s1 , s2 , . . . , sN ) and S is the maximal degree of the PES. The eigenfunctions of H0 noted Hermite functions

0 n (q)

are products of the one-dimensional normalized

ni (qi ):

0 n (q)

=

N Y

ni (qi ),

for all n = (n1 , . . . , nN ),

(2)

i=1

and the corresponding eigenvalues (in cm 1 ) are

En0

=

N X

1 (ni + )!i , 2 i=1

(3)

ni being the quantum number associated with the !i frequency. Let us define di the maximal degree of the Hermite function (see SI for the definition) associated with the i-coordinate in (2) and ⇧d the product space of size M of all quantum numbers defined by ⇧d =

N Y

[0, di ].

(4)

i=1

We seek the eigenfunctions

of H (1) as linear combinations of the eigenfunctions of H0

as (q) =

X

xn

0 n (q).

(5)

n2⇧d

By using the variational method on the eigenvalue problem for the operator defined in (1),

12

ACS Paragon Plus Environment

Page 13 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

we obtain the matrix form of the problem

Hf ull X = EX

(6)

where Hf ull is a M ⇥ M symmetric matrix and (E, X) is an eigenpair of Hf ull . An element of Hf ull (6) is calculated as

Hm,n = h

0 0 m |H n i

0 = Em

m,n

+h

0 0 m |V n i .

(7)

Solving Eq. (6) by iterative methods requires computation of fast matrix-vector products. Even if the Hamiltonian matrix is very sparse and the number of vavenumbers required (nv) is generally small, the size of the vectors, M , grows in a disjointed way. To overcome the curse of dimensionality, several approaches have been developped. These can be classified in two categories following how the basis is constructed: (i) fixed basis set (full VCI) ; (ii) adaptive basis set construction (variation-perturbation methods). The method used in this work concerns this second category. When an adaptive basis set construction is used, the eigenvalue computation of matrix H (6) is usually performed by a variational method combined with an iterative process. The matrix Hvp ⇢ Hf ull is diagonalized in embedded basis sets B (j) ⇢ B (j+1) ⇢ ⇧d ,

(8)

where B (j) is the basis set at iteration j. In the original variation-perturbation algorithm 62 configurations

0 m

(2) of B (j+1) are selected by imposing the following second order pertur-

bation criterion on the energy: h 0m |V 0n i2 > "VP , En0 0 En0 Em where

0 n

0 Em 6= 0,

(9)

2 B (j) and "VP is a given energetic threshold depending on the accuracy we want to 13

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 48

reach. The initial subspace B (0) includes the ground state and its excitations which represent the vibrational states to be determined. The convergence of the eigenvalue E (j) at iteration j is classically checked by evaluating |E (j)

E (j

1)

|

(10)

At the end of the variation-perturbation process Hvp is a Mvp ⇥ Mvp symmetric matrix (Mvp < M ). The energetic correction of all the (M

Mvp ) non-selected states 'n0 (q) ((2)) is later

estimated by perturbation theory. The sum B =

n 0 (q)

+ 'n0 (q) represents the global work

basis. Use of these approaches reduces the dimension of the variational problem dramatically 62,63 since the dimension of the active space to be diagonalized (Mvp ) is, now, under control. At each iteration only a small predetermined number (p) of vibrational states is needed to obtain the variational space containing the most important couplings. While this procedure converges the active space in a few (j) cycles and leads to matrices (Mvp ⇥ Mvp ) of only a few thousand terms for small systems,

Mvp = (p ⇥ j) + nv

(11)

this approach has yet to converge for systems larger than 20-25 atoms. For the sake of illustration, an example for a virtual system is reported in SI.

Despite the reduction in the size of the problem obtained with the use of the variationperturbation method, and before applying our approach to the sulfur polyaromatic compounds, it is important to finish this theoretical presentation by properly defining the dimension of the problem. As example, we have reported in Table 1 the values of the ⇧d space generated for the study of a 15 atom molecule and its dimer in terms of the maximal degree (d + N )! of the Hermite function d. The total number of elements of this binomial basis is N !d! 14

ACS Paragon Plus Environment

Page 15 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(see Table 1). It is very clear that the sizes of the matrices that will have to be generated in order to use the variational processes can very quickly reach several million lines and could not be used without the implementation of variation-perturbation methods for the study of medium size systems. Table 1: Size of the product space ⇧d (di =5, 6 or 7) for the study of a 15 atom molecule and its dimer. d

5

6

7

⇧d

1.086.000

8.145.060

53.524.680

⇧d

41.507.642

622.614.630

8.093.990.190

Reduction of the active vibrational space Here we show that a compromise is possible in order to further reduce the active space dimension that is required to converge the variational problem. Considering that all 3N

6

vibrational modes of a molecule are not uniformly coupled, two ‘truncated’ approaches were developed and tested in this work.

Truncation approach n 1 includes verifying the influence of the truncation of the basis B for studying a mode x (belonging to the T = 3N

6 ensemble) taking into account,

besides the mode being studied, only the X’ = x + a0 (a0 2 N) lower-frequency and the Y = [3N

6]

y higher-frequency normal modes (where y stands for the last non-selected

mode within the active space) (see Figure 2). Thus, one has a truncated problem dimension equal to T’ = X’ + Y < [3N

6]. A new ⇧0d product space of size M 0 is then defined where

⇧0d ⇢ ⇧d (For the sake of illustration, an example for a virtual system is reported in SI). The choice behind this truncation method relies on the fact that large amplitude modes, such as

15

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 48

torsion and rotation, are generally coupled to the elongation modes that act on lighter atoms, such as hydrogen. We are interested in these particular modes. Within this truncation, only p=100 vibrational states

0n 0 (q)

(2) by iteration j (11), strongly coupled to the state x for

which the vibrational characteristics are to be determined, are selected by a procedure that enriches the variational space. The energetic correction of all the (M 0

0 Mvp ) non-selected

states '0 n0 (q) is then estimated by perturbation theory. The sum B’ =

0n 0 (q)

+ '0 n0 (q) rep-

resents the new reduced global work basis.

Truncation approach n 2 was used to verify a hypothesis normally assumed to be true - that every x state is a priori more strongly coupled with its closest energetic neighbors. To do so, we have prepared an alternate selection criterion. Concerning the selection of the lower-frequency modes, we state a new criterion based on the selection of X”: in order to study a vibrational state x, one must select 2x” normal modes within the basis (x” being a positive integer) (see Figure 2). The selected 2x” normal modes within the truncated basis are exactly those that are energetically as close as possible to x. This basis has X”=1+2x” lower-frequency and Y higher-frequency states. The final dimension of this new basis is T”=X”+Y. A new ⇧00d product space of size M 00 is then defined where ⇧00d ⇢ ⇧d . Note that if x” equals zero, this is the same as studying the influence of high-energy stretching modes on x. This new procedure generates a new working basis denoted B” = ⌫1

y

x X’

(a) ⌫1

y

x

⌫3N

6

⌫3N

6

Y

X” (b)

+ '00 n0 (q).

Y

a0 x

00

00 n 0 (q)

x

00

Figure 2: Selection of normal modes which are stored in the expression of the PES : methodologies applied to realize the truncation of the vibrational basis a) truncation n 1 T’=X’+Y; b) truncation n 2 T”=X”+Y 16

ACS Paragon Plus Environment

Page 17 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Recalling Table 1 and with the goal of estimating the reduction of the size of the variational treatment, Table 2 reports, an example, the values of the reduced space (⇧0d or ⇧00d indi↵erently) generated for a 15 atom molecule and its dimer in the case (X 0 =7 and Y=21 i.e. 11 modes removed to the T full basis B) for the monomer and (X 0 =20 and Y=42 i.e. 22 modes removed to the T full basis B) for the dimer. It can be seen (Table 2) that the dimensions of the H matrix are greatly reduced with the use of these truncated bases but are still inaccessible for variation approaches.

Table 2: Size of the reduced product space ⇧0d or ⇧00d (di =5, 6 or 7) for the study of a 15 atom molecule and its dimer. d

5

6

7

⇧0d

237.336

1.344.904

6.724.520

⇧0d

9.657.648

109.453.344

1.078.897.248

Results and Discussion Monomers Benzothiophene (N=15) and dibenzothiophene (N=21) were chosen as references to test basis truncation. All the data concerning dibenzothiophene are reported in SI.

First, for each system, an accurate converged calculation of the variation-perturbation type was performed (T = 3N

6). The goal is to choose the subset of a ⇧d (4) giving the

best approximation of the eigenpairs. The first idea is to consider a complete active space Bd (di =5) including a maximum of a well-chosen 64 number of harmonic vibrational states

17

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0 n.

Second, the results obtained using the truncated bases (T 0 and T 00 ) are obtained. The

calculations were carried out using a variation-perturbation approach by fixing Mvp (max)= 25.000 vibrational states (2) and for nv=12 eigenvalues in all the calculations developed for the benzothiophene system. Third, the di↵erences obtained between the two procedures are calculated and are reported in Figure 3 (called D1 ) for truncation n 1) for the first and third modes (other modes among first twelve are provided in supporting information Figure S.5) and in Figure 4 (called D2 ) truncation n 2) for the fourth mode for benzothiophene (modes 5, 6 and 8 are provided in SI figure S.6). Of these twelve modes, modes 1, 2, 4, 5, 8, 10 and 12 occur in the Cs plane whereas modes 3, 6, 7 and 9 are out of the Cs plane. Since the outcomes of the test for the first nine modes of dibenzothiophene are comparable, see SI Figure S.7.

For truncation n 1, the main results are reported here in Figure 3. The information reported in these graphs represents the di↵erence obtained on the calculation of a target vibration frequency in the reduced basis B’ from the result that would be obtained using the same variation–perturbation method with the full basis B. For example, the analysis of the target mode 1 (X’=1, Figure 3) shows that the presence of only the higher-frequency !CH stretching modes suffices to describe the vibrational problem concerning this mode. This observation accords well with the hypothesis that only the large amplitude modes are strongly coupled with the elongation ones involving lighter atoms. The analysis of the other out-of-plane Cs deformation modes (2, 4, 5 and 8) leads to similar conclusions. The results of the calculation on mode 2 are closer to those obtained for mode 1, which seems logical since these two vibrational modes are de facto similar. Concerning modes 4 and 5, regardless of the number of selected states within the basis created from the X’ criterion, one can note the need to work with at least 13 modes on the Y basis in order to be as close as possible to the reference calculation. On the other hand, one observes that it is not necessary to take into account all the lower-frequency modes within the basis for the calculation to converge

18

ACS Paragon Plus Environment

Page 18 of 48

Page 19 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(D1 < 2cm 1 ). This will be confirmed shortly by the analysis of the second truncation. Another striking observation also emerges from this analysis: if it is desirable to introduce some of the wagging, scissoring and stretching modes in the basis, the calculation will converge satisfactorily only if all of these modes are introduced at the beginning. Concerning mode 8, it can be noted that no truncation seems capable of producing an acceptable convergence.

We have added the results for modes 10 and 12 (see SI Figure S.5) although they do not belong to the same low-frequency category considered in this study. These two modes lie at the interface between the low-frequency and the other modes. The main conclusion one can derive from them is that taking into account only the low frequency and !CH modes in the basis no longer suffices for convergence of the calculations (D1 < 2cm 1 ). These conclusions generally beyond mode 12.

The modes taking place in the plane of the molecule, i.e. modes 3, 6, 7, 9 and 11, again corroborate our hypothesis; particularly for low-frequency mode 3. In addition, the x0 modes that are close to mode 3 seem to have no e↵ect on the quality of the converged result. The descriptions of modes 6 and 7 show again that the results are not very dependent on the basis truncation unless all the !CH modes were previously selected. Analyzing the results for modes 9 and 11 leads to the same conclusion, except that a strong truncation of Y would cause a much more important deviation of the results from those observed for modes 3, 6 and 7.

The analysis of the results concerning the second truncation confirms even more the observations just presented: only the results obtained for mode 4 are reported here (modes 5, 6 and 8 are provided in SI - Figure S.6). Representations of these results were added as a second layer to those obtained beforehand (Figure 4); the two situations are connected by a series of points they share in common, i.e. when a0 =x” or a0 =x

19

ACS Paragon Plus Environment

1 (see SI Figure S.3). The

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

main information we can derive from these studies is that no low-frequency mode requires its closest neighbors to ensure a good convergence of the results. Every argument indicates that these modes are only coupled to those involving !CH stretching. Only modes 4 and 5 seem to slightly depend on each other.

In summary, variational-perturbational studies using a truncated basis taking into account the closest neighbors of the low-frequency modes seems not to be an essential requirement for describing their vibrational properties as opposed to their coupling with the highfrequency stretching !CH modes. Finally, variational-pertubational studies using a truncated basis (truncation n 1) were performed for systems 1a (benzothiophene), 1b (dibenzothiophene), 1c (4,6-dimethyldibenzothiophene) and 1d (4-methyldibenzothiophene). The same was done for their dimers and compared against experimental results obtained by IR and Raman spectroscopy.

(a)

20

ACS Paragon Plus Environment

Page 20 of 48

Page 21 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(b)

Figure 3: Results obtained using truncation n 1 for (a) first (in Cs plane) and (b) third (out of Cs plane) mode of a benzothiophene molecule. In the z axis reports the values of the di↵erence (D1 ) obtained on calculating a frequency without any truncation and the same frequency calculated in base B’ (T’=X’+Y; X’ axis x; Y axis y; in accordance with the notation used in Figure 2).

PA infrared and Raman spectra of 1b, 1c and 1d were interpreted in light of the quantum chemistry calculations. Data for 1a were taken from the literature. The low-wavenumber and “fingerprint” regions (comprising, together, the 20–2000 cm

1

range) were studied in de-

tail in these investigations. In addition to the predicted infrared-active bands, the PA spectra of these hydrocarbons display many features due to overtones and combinations. Moreover, numerous Raman-active (gerade) vibrations also give rise to bands in the PA spectra, which thus convey considerable information regarding the structures of these compounds.

21

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4: Result using the truncation n 2 for fourth vibration mode of a benzothiophene molecule. In the z axis reports the values of the di↵erence (D2 ) obtained on calculating a frequency without any truncation and the same frequency calculated in base B” (T”=X”+Y; X” axis x; Y axis y; in accordance with the notation used in Figure 2). Low wavenumber region ( 100–700 cm 1 ): both Raman and infrared spectroscopy results for 1a and 1b are reported in Tables 3 and 4. Raman and PA infrared results for 1c and 1d are given in the SI (Tables S.7 and S.9).

For the benzothiophene molecule, the reported Raman spectrum has no significant features in this region. The IR spectrum is dominated by three intense bands at 471, 560 and 690 cm

1

which are characteristic of wagging modes. Two other combination modes are

IR-active. They are not easily detected since they have very low intensity and also due to the fact that the ⌫2 + ⌫5 mode has a signature very close to that of the active fundamental mode ⌫9 . The analysis of the active modes of the dibenzothiophene molecule (Figure 5) 22

ACS Paragon Plus Environment

Page 22 of 48

Page 23 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

shows that two of the three most intense wagging modes (obs 497 and 704 cm 1 ), which are, incidentally, characteristic of benzothiophene, can be clearly distinguished and do not allow one to easily separate the molecules from each other. However, one should note that a characteristic shift towards higher frequencies is found for these two modes of the dibenzothiophene system. Also note, in comparison to benzothiophene, an extinction (near 560 cm 1 ) of the wagging signature because of symmetry induction and a new intense characteristic vibration located around 615 cm 1 . Moreover, the IR spectrum of this molecule for the region between 100 and 450 cm

1

is more complex and richer in features. As a matter of

fact, other stronger and medium intensity modes signify the presence of the second benzyl moiety within the dibenzothiophene molecules. All other active modes have very low intensities, which again agrees with our calculations. The analysis of the eigenfunctions generated by the variational-perturbational process shows that these modes are mainly developed on the basis of vibrational states that are characteristic of torsion (⌫3 (a1 ) 213 cm 1 ; calc 217 cm 1 ), breathing (⌫4 (b1 ) 228 cm 1 ; calc 221 cm 1 ) and (⌫6 (a1 ) 404 cm 1 ; calc 416 cm 1 ) and wagging (⌫8 (b1 ) 425 cm 1 ; calc 433 cm 1 ) combination modes of the benzenic rings.

The Raman signature of the dibenzothiophene molecule is also more complex than that of benzothiophene. In the present context, this fact introduces only limited information since most of the active Raman modes are already the most intense and characteristic of the IR spectrum. All the other active modes have very low intensities which are again in agreement with our calculations. It is worth noting the confirmation provided by the Raman results concerning the shoulder observed in the IR spectrum around 419 cm 1 . Similarly to the benzothiophene molecule two combination modes are foreseen by the calculations that can be, mainly for the ⌫2 + ⌫8 mode, reasonably assigned experimentally.

In addition to some frequency shifts, the spectra of both 4-methyldibenzothiophene and 4,6-dimethyldibenzothiophene exhibit important changes in the intensities of the vibrations

23

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

in the 100 - 500 cm

1

Page 24 of 48

region as compared with the spectrum of dibenzothiophene. The

clearest change occurs in the 400-460 cm

1

region: the presence of substituents induces a

significant change in intensity, ranging from very intense modes down to nearly zero (notably for the ⌫10 mode of 4-methyldibenzothiophene). These di↵erences are even more pronounced in the Raman spectra. Whilst the presence of a single methyl group does not significantly change the Raman signature of the dibenzothiophene molecule, a second substitution in the 4,6 position does so very pronouncedly. The main change concerns the most characteristic mode of the di-substitution found at 598 cm 1 , which represents an in-phase breathing mode of the three rings. Table 3: Low wavenumber Raman and infrared spectra of benzothiophene (Cs symmetry ) Observedb Raman Infrared 192 vw 340 vw 387 vw

491 vw 526 vw

192 vw 198-212 w 340 vw 405-411 vw 471 m 491-497 vw 558-564 m 669 vw 690 s

Calculation ! a,c

Mode ⌫1 (A”) ⌫2 (A”) ⌫3 (A’) ⌫1 + ⌫2 (A’) ⌫4 (A”) ⌫5 (A”) ⌫6 (A’) ⌫7 (A’) ⌫8 (A”) ⌫9 (A’) ⌫2 + ⌫5 (A’) ⌫10 (A”)

196 (0.22) [0.27] 204 (4.11) [0.12] 351 (2.82) [2.11] 427 (3.45) [0.77] 482 (5.88) [0.53] 500 (0.30) [10.64] 538 (0.25) [6.54] 578 (5.88) [0.70] 683 (2.32) [1.47] 710 (34.70) [0.23]

a

⌫d 194 (0.20) 199 (2.97) 345 (2.08) 391 (0.03) 414 (2.50) 478 (6.02) 494 (0.29) 533 (0.36) 563 (5.24) 675 (1.70) 679 (0.10) 692 (28.22)

Infrared intensities (in parentheses) in km/mol and Raman intensities [in brackets] in b A4 /AMU. Experimental data from ref. 65 Relative intensities: vw = very weak, w c = weak, m = medium, s = strong, vs = very strong, sh = shoulder. Unscaled hard monic frequencies calculated using !B97X-D/6-311++G**. Anharmonic frequencies !B97X-D/6-311++G** calculated using the P Anhar algorithm. 60

24

ACS Paragon Plus Environment

Page 25 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 4: Low wavenumber Raman and infrared spectra of dibenzothiophene (C2v symmetry ) Observedb Raman Infrared

Mode

138 w 171 w 217 m 232 w

⌫1 ⌫2 ⌫3 ⌫4

282 w 287 sh 410 s 419 m 437 vw 496 w 504 vw

135 m 213 m 228 m 269 vw 282 vw 404 s 419 sh 425 s 436 w 460 vw 497 s 512 w 565 w 615 m 667 w

703 s

704 vs

Calculation ! a,c

(b1 ) (a2 ) (a1 ) (b1 )

102 135 217 221

(1.54) (0.00) (0.85) (1.46)

[0.28] [0.07] [1.25] [1.03]

⌫d 101 133 208 223

(0.66) (0.01) (1.23) (1.42)

⌫5 (a2 )

279 (0.00) [2.90]

274 (0.00)

⌫6 ⌫7 ⌫8 ⌫9

416 (2.83) [13.05] 428 (0.14) [2.72] 433 (7.15) [0.01] 444 (0.00) [0.06]

409 425 425 438

(2.80) (0.15) (5.06) (0.00)

505 (0.28) [4.76] 511 (1.23) [4.67] 510 (2.74) [0.03]

499 497 504 562 568 627 666 675 713

(0.22) (1.48) (1.18) (0.12) (0.00) (4.62) (0.15) (0.01) (6.43)

(a1 ) (b2 ) (b1 ) (a2 )

⌫10 (1 ) ⌫12 (b2 ) ⌫11 (b1 ) ⌫2 + ⌫8 (b2 ) ⌫13 (a2 ) ⌫14 (b2 ) ⌫17 (b1 ) ⌫3 + ⌫10 (b2 ) ⌫15 (b2 )

580 (0.00) [0.24] 635 (5.09) [0.04] 730 (0.82) [0.00] 724 (6.53) [5.58]

a

Infrared intensities (in parentheses) in km/mol and Raman intensities [in brackets] in A4 /AMU . b Experimental data from this work and ref. 30 Relative intensities: vw = very weak, w c = weak, m = medium, s = strong, vs = very strong, sh = shoulder. Unscaled hard monic frequencies calculated using !B97X-D/6-311++G**. Anharmonic frequencies !B97X-D/6-311++G** calculated using the P Anhar algorithm. 60

25

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 48

Figure 5: Far-infrared PA spectra of dibenzothiophene

Dimers Analysis of the interaction energies components Characterization of the spectral signatures of “isolated” thiophene molecules constitutes the first part of this study. The second aim of this work exists in the identification and characterization of the IR and Raman signatures of the vibration modes at very low wavenumbers which are characteristics of associated systems such as dimers. In order to foresee and identify the spectral signatures typical of these interactions, we have first identified all the stable forms of benzothiophene (and its derivatives) dimers, i.e., 2a d -benzothiophene, 2b d dibenzothiophene, 2c d -4-methyldibenzothiophene and 2d d -4,6-dimethylbenzothiophene. Among all the possible dimer configurations, eight types of benzothiophene dimers have been calculated; these are represented in Figure 6. In all tables and figures, the letter d means dimer. The geometrical parameters are summarized in Table 5. Our interest lies in ⇡-stacking interactions, hence we have considered the molecules placed such that ⇡ stacking interactions are favored. 66,67 26

ACS Paragon Plus Environment



Page 27 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Although ⇡

⇡ stacking interactions have been proven to be predominant, electrostatic

interactions also induce specific configurations. 18 In order to validate the calculation method used, additional calculations (geometrical parameters and interaction energies) using SAPTDFT were also performed. The energies for the di↵erent configurations of the benzothiophene molecule are reported in Table 5. The energies for all the configuration for the three other systems (2b-2c), and the graphical representations and optimized structures can be found in the SI (Pages S29-S32).

Figure 6: Configurations of benzothiophene dimers optimized at the !B97X-D/6-311++G** level of theory. Sulphur atoms are shown in yellow.

Figure 7: Structural parameters of benzothiophene dimers

d = Intermolecular distance between thiophene rings center of benzothiophene dimers. ↵ = Angle between the molecules inside x, z plane. = Angle between the molecules inside y, z plane. = Angle of rotation in z axis.

27

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 48

Eight conformations were highlighted for the benzothiophene dimer species. Benchmark energy calculations were performed using the DFT-SAPT(PBE0)/aug-cc-pvtz level of theory. The !B97X-D/ 6-311++G** calculations, from which the vibrational calculations were performed, lead to the same conclusions concerning the stability of these eight conformers (average deviation on the order of 1.25 kcal/mol). Regardless of the theory level that is employed, the two most stable conformations are D and E. In these conformations, the sulfur atoms are opposed to each other. This is also the case for the majority of the other systems studied in this work.

Before commenting on the results obtained by these contributions, one should briefly recall the basis of the SAPT method used here. The van der Waals interactions, responsible for the cohesion of the dimers reported in this work, are clearly distinct to other types of long-range interactions. In the intermediate region that includes the minimum, it is very difficult to unequivocally decompose the total interaction energy into di↵erent contributions. Several decomposition schemes based on the Rayleigh-Schr¨odinger perturbation theory have been proposed over the years. The principle of these approaches is the separation of the di↵erent contributions, such as electrostatic, induction and dispersion. On the other hand, the greatest drawback is the difficulty regarding the convergence of the perturbative series (mainly the polarisation part, i.e., induction + dispersion), due to their truncation at low orders or their incompleteness. SAPT has been proposed to address this. 36 In this approach, each of the aforementioned components should also include one exchange correction in order to take the antisymmetric nature of the wave function of the created complex into account. A brief review of the SAPT method and the main contributions / notations that are used in Table 6 are reported in SI (Fig S.12).

28

ACS Paragon Plus Environment

Page 29 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 5: (a) Interaction energies of benzothiophene dimers (kcal/mol) in 8 different configurations using !B97X-D/6-311++G**, and compared with other methods. (b) Intermolecular distance between thiophene rings center of benzothiophene dimers (˚ A) (c) Angle between ring planes of benzothiophene dimers ( ) according to Figure 7 and calculated at the !B97X-D/ 6-311++G** level of theory. Energies taking into account zero-point energy (ZPE) and Basis Set Superposition Error (BSSE) corrections. (a) A

B

C

D

E

F

G

H

-6.41 -5.33 -5.46

-7.05 -5.73 -5.80

-6.83 -5.56 -5.72

-7.11 -5.95 -5.68

-7.23 -5.83 -5.13

-7.01 -5.87 -5.40

-7.02 -5.72 -5.10

-7.00 -5.63 -4.76

A

B

C

D

E

F

G

H

3.94

3.81

5.28

3.72

4.49

3.56

3.74

4.37

DFT-SAPT(PBE0)/aug-cc pVTZ

3.98

3.80

5.16

3.75

4.53

3.65

3.70

4.30

B971/6-31+G(d,p)-DCPa

3.99

3.84

4.65

3.83

5.03

3.66

3.75

3.70

A

B

C

D

E

F

G

H

2.94 3.59

1.73 4.23

1.06 3.14

0.44 0.51

2.23 5.07 68.92

4.38 13.28 100.06

2.05 9.03 51.71

6.37 5.64 77.25

!B97x-D/6-311++G** DFT-SAPT(PBE0)/aug-cc pVTZ B971/6-31+G(d,p)-DCPa

(b) !B97x-D/6-311++G**

(c) ↵

a

ref 25

One should note that the SAPT approach used for the study of non-polar complex systems (for which van der Waals forces are predominant) brings up little information and it is unnecessary to develop the perturbation beyond the second order. The so-obtained results (truncated at the order 2) for the four dimers B, D, E and F of the benzothiophene system are reported in Table 6. The very small numerical di↵erences observed show that the nature 29

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 48

of the interaction energies for di↵erent conformers does not make it possible to prove the existence of one as opposed to the others. Recalling the Introduction to this article, one could accordingly propose the concomitant presence of all these conformations as representative of inhomogeneous mixtures of asphaltenes, justifying the analysis made herein. Table 6: Interaction energies components in kcal/mol of benzothiophene dimers using DFT-SAPT(PBE0)/aug-cc pVTZ

Ecorr inter Ecorr intra Eint

B

D

E

F

-12.00 -0.93 -5.73

-11.78 -0.99 -5.95

-12.00 -0.84 -5.83

-10.57 -0.92 -5.87

In order to complete the analysis of the interaction energies, the Non-Covalent Interactions (NCI) formalism developed by Johnson and coworkers 68 was also used to analyze the energetic di↵erences between all stable dimer forms. This formalism is based on the reduced electronic density gradient (s), a fundamental quantity in DFT theory which describes the deviation of the electronic density of a system (⇢) from the homogeneous case. Plotting s versus ⇢ reveals the basics of intramolecular interactions: for small ⇢ and large s, one can identify the exponential decay of the density far from the nuclei. For small s and large ⇢, the points correspond to the covalent bond between the atoms. Situations where ⇢ and s are both small are diagnostic for non-covalent interactions. To distinguish between attractive and repulsive non-covalent interactions, one needs to consider the Laplacian of the density (r2 ⇢).

The components of the decomposition along the three coordinates are the three eigenvalues

i

of ⇢’s Hessian matrix in such a way that r2 ⇢ =

Near the nuclei, there is a maximum for ⇢ and all between bonded atoms, the situation

3

> 0 and

i

1

+

2

+

3

is true for

3

2

1.

are negative. For the interatomic regions 1

and

2

< 0 exists. In this formalism,

for covalent interactions, the negative contributions are dominant and the Laplacian is then 30

ACS Paragon Plus Environment

Page 31 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

negative. On the other hand, for non-covalent interactions, the Laplacian in the interatomic region is positive, irrespective of the bonding or non-bonding character of the interaction. Bonding interactions are, in this way, identified by a negative value of

2

whereas a positive

value indicates non-bonding contacts. In other words, the sign of this eigenvalue is used to distinguish bonded (

2

< 0) and non-bonded (

2

> 0) interactions. ⇢ itself indicates the

strength of the interaction, as can be seen for the two most stable conformers D and E from Figures 8a and 8b respectively. Pictures concerning all the other conformers are reported in SI.

(a)

(b)

Figure 8: Plots of the electronic reduced gradient density against the density ⇢ times the sign of 2 - Cases of the conformers D (a) and E (b). In Figure 8, the low density, high gradient regions correspond to the exponential decay of the electronic density away from the nuclei of the molecules. The region of gradient values around 0.5 arbitrary units and density values around ± 0.2 au correspond to the covalent interactions: in the positive region this is related to the repulsive character of the cova31

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 48

lent intramolecular bonds and the negative region is correlated with attractive character. Moreover, the plots of these surfaces allow one to visualize the regions of space where the interactions are stronger. The color mapping indicates their repulsive character (towards red) or attractive (towards green). The intra-ring nuclear repulsion can be clearly distinguished from the interaction surfaces between the molecules of the dimer system in this way.

From these plots, one can deduce that the conformations of the dimers have, to a limited extent, some influence on the non-covalent interactions. Although very little can be said from the reduced density gradient versus the density plots, the non-covalent surfaces show singular features for the studied conformers: conformer F has attractive interactions between the hydrogen atoms of one molecule and the carbon atoms of the other while conformer H has regions where no interaction takes place or where they mutually cancel (indicated by the holes on the surface). Moreover, analysis of the positioning of the intra-ring nuclear repulsion interactions shows that systems C, D and E have the surfaces of each benzyl or thiophene ring in opposition to each other. This increased interaction between the planes might indicate that the translational (sliding) vibration modes of one plane over the other should be sliding shifted to higher wavenumbers.

Interpretation of Spectra Very Low wavenumber region ( 30–100 cm 1 ) : The first conclusion to be drawn from our calculations for all conformations for all species is that the very low wavenumber features in the spectra for all of the dimers are very similar to each other.

This observation is also consistent with our spectra recorded between 5 and 26 cm

1

for

the 1b, 1c and 1d systems (see Figure 9). The assigment of the six libration modes of the most stable conformers for the four species (2a-2d ) is reported in Table 7. The parame32

ACS Paragon Plus Environment

Page 33 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ters of the variation-perturbation calculations used are reported in Table 8 (truncation n 1). The analysis of the vibrational wavefunctions obtained from the diagonalization of the Hvp matrix show that the modes with the lowest wavenumbers, i.e., the sliding mode along the xy plane (Figure 7) and the rotation around the main axis perpendicular to the plane of the isolated molecules (y direction) are mainly coupled; these modes are also strongly coupled with several CH stretching modes. On the other hand, the other libration modes and, particularly, the sixth corresponding to the ↵-type rocking mode are coupled, at the same time, to the stretching and wagging !CH modes. In particular, the results concerning the dibenzothiophene molecule are reported in Table 9 and compared to the only experimental data available in the literature and also to our PA measurements. All the results that were gathered from the literature do not lead to the same frequency, partly due to the di↵erent methodologies in use, but rather to a dispersion (around the frequency found by us). This is unsurprising since the measurement conditions di↵er among these experiments. The calculations performed for the dibenzothiophene dimer are in good agreement with the ensemble of experimental observations. They confirm that the very low-intensity bands observed experimentally are indeed the signatures of the intermolecular modes. The results also show that the rocking modes, found between 60 and 100 cm 1 , are not strongly dependent on the species under study, and have more intense signatures than the libration modes. The rotation and sliding modes are much less intense and their locations are more dependent on the type of molecule. They are also more difficult to characterize in all types of spectroscopy: however, we have recently identified weak infrared absorption bands near 20 cm

1

that may

correspond to the sliding mode in spectra acquired using coherent synchrotron radiation (Figure 9). Finally, the presence of one or several substituents at di↵erent positions on the molecule can also change the intensity of the signatures of the libration modes significantly, particularly for those implying, through vibration, the movement of the substituents towards each other, which gives rise to new intermolecular interactions.

33

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 48

Previously, we showed the extent to which these interactions play a significant role in the aggregation of molecules constituting the asphaltenes and the importance of characterizing them. 67 Among these interactions, a non-negligible part consists of intermolecular interactions. Characterizing them through spectroscopy is then crucial for understanding aggregation phenomena of asphaltene mixtures and, thus, estimating the energies required to separate them during the refining process and avoid the asphaltene sedimentation phenomenon widely observed in the pipelines. Our work using molecular dynamics sheds some light on the nature of the intermolecular interactions that can take place in the mixture for di↵erent heteroatoms, di↵erent lateral chains, etc. Among our future aims is the identification of characteristic features in vibrational spectra which permit identification of each interacting species for di↵erent heteroatoms and lateral chains that may be present in these aggregated systems.

Table 7: calculated intermolecular vibrations modes of the four most stable dimers 2a-2d

Benzothiophene ! a (I)b

Mode ⌫1 ⌫2 ⌫3 ⌫4 ⌫5 ⌫6 a

: cm

17 18 45 67 74 91 1

(0.02) (0.08) (0.06) (0.32) (0.02) (0.06) b

⌫a 17 18 39 58 67 88

Dibenzothiophene ! a (I)b 18 28 36 70 89 93

(0.05) (0.01) (0.09) (0.32) (0.01) (0.56)

⌫a 18 24 35 66 82 92

4-methyl dibenzothiophene ! a (I)b ⌫a 22 36 37 65 81 96

(0.01) (0.17) (0.01) (0.07) (0.03) (0.06)

: km/mol

34

ACS Paragon Plus Environment

22 32 31 55 67 85

4,6-dimethyl dibenzothiophene ! a (I)b ⌫a 17 22 41 74 83 84

(0.00) (0.01) (0.01) (0.52) (0.07) (0.01)

17 22 38 66 69 78

Page 35 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 8: Parameters used in the variation-perturbation process (values of X’ and Y as defined in Figure 2 parameters and dimension of the variational problem Mvp ); p=100; nv=12 (see Eq. (11)) for each intermolecular mode of the four most stable dimers 2a-2d

Benzothiophene Mode

type

X’ ; Y ; Mvp

⌫1 ⌫2 ⌫3 ⌫4 ⌫5 ⌫6

x,y

1 ; 7 ; 2500 3 ; 7 ; 2500 5 ; 7 ; 3000 7 ; 12 ; 3600 7 ; 12 ; 3500 9 ; 12 ; 4200

x,y ↵ ↵

Dibenzothiophene type

X’ ; Y ; Mvp

x,y x,y ↵ ↵+!

1 3 5 7 7 11

; ; ; ; ; ;

10 10 10 15 15 15

; ; ; ; ; ;

3200 4100 4500 8800 9300 15000

4-methyl dibenzothiophene type X’ ; Y ; Mvp

x,y x,y ↵ ↵+!

1 ; 14 ; 8100 3 ; 14 ; 9000 5 ; 14 ; 10500 7 ; 20 ; 12700 7 ; 20 ; 12200 11 ; 20 ; 17200

4,6-dimethyl dibenzothiophene type X’ ; Y ; Mvp

x,y x,y ↵ ↵+!

1 ; 17 ; 10300 3 ; 17 ; 15100 5 ; 17 ; 16600 7 ; 23 ; 20500 7 ; 23 ; 17400 11 ; 23 ; 24900

Table 9: Experimental and calculated intermolecular vibration modes (cm 1 ) of dibenzothiophene dimer (B) IRa

!B97x-D/6-311++G** Mode ⌫1 ⌫2 ⌫3 ⌫4 ⌫5 ⌫6 a

! (I)c [IRaman ]d



18 28 36 70 89 93

18 24 35 66 82 92

(0.05) (0.01) (0.09) (0.32) (0.01) (0.56)

[0.62] [7.96] [9.15] [0.51] [0.51] [3.33]

: liquid/crystal ref. 30

b

Ramana

IRb

⌫ ⇠ 20 vw 53 w 72 w 89 vw 104 - 101 w 109 w

: this work (crystal)

35

c

48 - 49 63 79 - 82 103 - 106 126

: km/mol

ACS Paragon Plus Environment

d

56 w 69 vw 78 vw 88 vw

: A4 /AMU

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a) dibenzothiophene

(b) 4,6-dimethyldibenzothiophene

(c) 4-methyldibenzothiophene

Figure 9: Far-infrared spectra of dibenzothiophene 1b, 4,6-dimethyldibenzothiophene 1c and 4-methyldibenzothiophene 1d . Data were acquired using coherent synchrotron radiation and a Bruker IFS 125HR FT-IR spectrometer.

36

ACS Paragon Plus Environment

Page 36 of 48

Page 37 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Conclusion In summary, the present work shows that: 1) solving the vibrational Schr¨odinger equation for isolated molecules allows one to identify, for molecules containing the sulfur heteroatom, the vibration modes for di↵erent molecular moieties as well as their singular spectroscopic properties; 2) utilizing both IR and Raman spectroscopy is very useful towards the full identification of these modes; 3) variation-pertubation calculations in the 100-700 cm

1

re-

gion make it possible to complete the assignment of the modes of each monomer (which was previously initiated in the literature); 4) in the same region, vibrational modes of the dimers correspond to those of the isolated molecules making it possible to determine, based solely on these modes, whether the molecules are interacting or not; 5) these intermolecular interaction modes are expected only at very low wavenumbers, and 6) in the region below 100 cm 1 , each molecule containing di↵erent configurations has characteristic signatures in the vibrational spectrum: these signatures arise from the sliding or rotation modes of one molecular plane with respect to another and can be in- or out-of-phase. Even though no general rule could be established concerning the behavior of these modes for di↵erent species or configurations, it is clear that their possible existence can serve as a marker for the presence or absence of intermolecular aggregation.

Additionally use of the NCI formalism on a wide range of molecules cointaining di↵erent substituents and di↵erent heteroatoms should enable us to characterize the nature of these interactions in the future, and to quantify the energies involved in these interactions. This will make it possible to predict the behaviors for di↵erent molecular architectures and heteroatom arrangements.

Finally, the description of the intermolecular modes chosen in this work is a first step of a project aimed at fully characterizing the aggregation interactions. In future work, we plan to a) perform solid state calculations, which are closer to actual experimental conditions 37

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 48

since these interaction modes are linked to the deformations of the crystalline network of the solid (related results were already published by our group concerning the vibrational modes near 100 cm

1

for a family of acenes 29 ); and b) study families of molecules with di↵erent

chemical substituents and di↵erent heteroatoms (N, O and S).

Acknowledgement Part of the work described in this paper was performed at the Canadian Light Source, which is supported by the Natural Sciences and Engineering Research Council of Canada, the National Research Council of Canada, the Canadian Institutes of Health Research, the Province of Saskatchewan, Western Economic Diversification Canada, and the University of Saskatchewan. Computer time for this study was provided by computing facilities MCIA (M´esocentre de Calcul Intensif Aquitain) at the Universit´e de Bordeaux and at the Universit´e de Pau et des Pays de l’Adour.

Support from the sponsors of the NSERC IRC in Petroleum Thermodynamics (the Natural Sciences and Engineering Research Council of Canada (NSERC), Alberta Innovates Energy and Environment Solutions, BP Canada, ConocoPhillips Canada Resources Corp., Nexen Energy ULC, Shell Canada Ltd., Total E&P Canada Ltd., and the Virtual Materials Group) is also gratefully acknowledged. We thank Dr C. Darrigan and Dr G. Vallverdu for their constructive suggestions on graphical techniques.

Supporting Information Available The following files are available free of charge. The Supporting Information is available free of charge on the ACS Publication website at DOI : insert DOI number Variation-Perturbation principe; The two truncations are connected by a series of points

they share in common; The benzothiophene and dibenzothiophene molecules; The 4-methyldibenzothiophen 38

ACS Paragon Plus Environment

Page 39 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

and 4,6-methyldibenzothiophene molecules; The dimers molecules; Brief informations onto the symmetry-adapted perturbation theory; Non-Covalent Interactions (NCI) results; Hermite functions.

References (1) Waldo, G. S.; Carlson, R.; Moldowan, J.; Peters, K. E.; Penner-hahn, J. E. Sulfur Speciation in Heavy Petroleums: Information from X-ray Absorption Near-Edge Structure. Geochimica et Cosmochimica Acta 1991, 55, 801–814. (2) Mitra-Kirtley, S.; Mullins, O. C.; Ralston, C. Y.; Sellis, D.; Pareis, C. Determination of Sulfur Species in Asphaltene, Resin, and Oil Fractions of Crude Oils. Applied spectroscopy 1998, 52, 1522–1525. (3) Lobodin, V. V.; Robbins, W. K.; Lu, J.; Rodgers, R. P. Separation and Characterization of Reactive and Non-Reactive Sulfur in Petroleum and Its Fractions. Energy & Fuels 2015, 29, 6177–6186. (4) Katz, H. E. Recent Advances in Semiconductor Performance and Printing processes for Organic Transistor-Based Electronics. Chemistry of Materials 2004, 16, 4748–4756. (5) Ong, B. S.; Wu, Y.; Li, Y.; Liu, P.; Pan, H. Thiophene Polymer Semiconductors for Organic Thin-Film Transistors. Chemistry–A European Journal 2008, 14, 4766–4778. (6) Vivas, M. G.; Nogueira, S. L.; Silva, H. S.; Barbosa Neto, N. M.; Marletta, A.; SereinSpirau, F.; Lois, S.; Jarrosson, T.; De Boni, L.; Silva, R. A. Linear and Nonlinear Optical Properties of the Thiophene/Phenylene-Based Oligomer and Polymer. The Journal of Physical Chemistry B 2011, 115, 12687–12693. (7) Silva, H. S.; Nogueira, S. L.; Manzoli, J. E.; Barbosa Neto, N. M.; Marletta, A.; SereinSpirau, F.; Lere-Porte, J. P.; Lois, S.; Silva, R. A. Controlling Bandgap Energy and Mul39

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tivibronic Modes of a Poly (2,5-thiophene-1,4-dialkoxyphenylene)Derivative by Gamma Photons. The Journal of Physical Chemistry A 2011, 115, 8288–8294. (8) Kim, F.; Ren, G.; Jenekhe, S. A. One-dimensional Nanostructures of ⇡-Conjugated Molecular Systems: Assembly,Pproperties, and Applications from Photovoltaics, Sensors, and Nanophotonics to Nanoelectronics. Chemistry of Materials 2010, 23, 682–732. (9) Kishita, A.; Takahashi, S.; Kamimura, H.; Miki, M.; Moriya, T.; Enomoto, H. Upgrading of Bitumen by Hydrothermal Visbreaking in Supercritical Water with Alkali. Journal of the Japan Petroleum Institute 2003, 46, 215–221. (10) Samokhvalov, A. Heterogeneous Photocatalytic Reactions of Sulfur Aromatic Compounds. ChemPhysChem 2011, 12, 2870–2885. (11) Chen, X.; Gao, J.; Lu, Y.; Meng, H.; Li, C. Acylation Desulfurization of Heavy Cracking Oil as a Supplementary Oil Upgrading Pathway. Fuel Processing Technology 2015, 130, 7–11. (12) Strausz, O. P.; Mojelsky, T. W.; Lown, E. M. The Molecular Structure of Asphaltene: an Unfolding Story. Fuel 1992, 71, 1355–1363. (13) Liu, P.; Shi, Q.; Chung, K. H.; Zhang, Y.; Pan, N.; Zhao, S.; Xu, C. Molecular Characterization of Sulfur Compounds in Venezuela Crude Oil and its SARA Fractions by Electrospray Ionization Fourier Transform Ion Cyclotron Resonance Mass Spectrometry. Energy & Fuels 2010, 24, 5089–5096. (14) Karimi, A.; Qian, K.; Olmstead, W. N.; Freund, H.; Yung, C.; Gray, M. R. Quantitative Evidence for Bridged Structures in Asphaltenes by Thin Film Pyrolysis. Energy & Fuels 2011, 25, 3581–3589. (15) Coelho, R. R.; Hovell, I.; Rajagopal, K. Elucidation of the Functional Sulphur Chemi-

40

ACS Paragon Plus Environment

Page 40 of 48

Page 41 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

cal Structure in Asphaltenes using First Principles and Deconvolution of Mid-Infrared Vibrational Spectra. Fuel processing technology 2012, 97, 85–92. (16) AlHumaidan, F. S.; Hauser, A.; Rana, M. S.; Lababidi, H. M. Impact of Thermal Treatment on AsphalteneFunctional Groups. Energy & Fuels 2016, 30, 2892–2903. (17) Widany, J.; Daminelli, G.; Di Carlo, A.; Lugli, P.; Jungnickel, G.; Elstner, M.; Frauenheim, T. Electronic Band Structure and Intermolecular Interaction in Substituted Thiophene Polymorphs. Physical Review B 2001, 63, 233204. (18) Liu, H.; Br´emond, E.; Prlj, A.; Gonthier, J. F.; Corminboeuf, C. Adjusting the Local Arrangement of ⇡-Stacked Oligothiophenes through Hydrogen Bonds: a Viable Route to Promote Charge Transfer. The journal of physical chemistry letters 2014, 5, 2320– 2324. (19) Speight, J. G. The Chemistry and Technology of Petroleum, 5th ed.; CRC press, Taylor & Francis Group, 2014. (20) Mullins, O. C.; Sabbah, H.; Eyssautier, J.; Pomerantz, A. E.; Barr´e, L.; Andrews, A. B.; Ruiz-Morales, Y.; Mostowfi, F.; McFarlane, R.; Goual, L. et al. Advances in Asphaltene Science and the Yen–Mullins Model. Energy & Fuels 2012, 26, 3986–4003. (21) Sheremata, J. M.; Gray, M. R.; Dettman, H. D.; McCa↵rey, W. C. Quantitative Molecular Representation and Sequential Optimization of Athabasca Asphaltenes. Energy & Fuels 2004, 18, 1377–1384. (22) Silva, H. S.; Sodero, A. C.; Korb, J.-P.; Alfarra, A.; Giusti, P.; Vallverdu, G.; B´egu´e, D.; Baraille, I.; Bouyssiere, B. The Role of Metalloporphyrins on the Physical-Chemical Properties of Petroleum Fluids. Fuel 2017, 188, 374–381. (23) Speight, J. G. The Desulfurization of Heavy Oils and Residua, 2nd ed.; CRC Press, 1999. 41

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(24) Groenzin, H.; Mullins, O. C. Asphaltene Molecular Size and Structure. The Journal of Physical Chemistry A 1999, 103, 11237–11245. (25) Mackie, I. D.; DiLabio, G. A. Importance of the Inclusion of Dispersion in the Modeling of Asphaltene Dimers. Energy & Fuels 2010, 24, 6468–6475. (26) Morgan, R. S.; Tatsch, C. E.; Gushard, R. H.; Mcadon, J. M.; Warme, P. K. Chains of Alternating Sulfur and ⇡-Bonded Atoms in Eight Small Proteins. International journal of peptide and protein research 1978, 11, 209–217. (27) Tauer, T. P.; Derrick, M. E.; Sherrill, C. D. Estimates of the ab initio Limit for Sulfur-⇡ Interactions: the H2 S-Benzene Dimer. The Journal of Physical Chemistry A 2005, 109, 191–196. (28) Tsuzuki, S.; Honda, K.; Azumi, R. Model Chemistry Calculations of Thiophene Dimer Interactions: Origin of ⇡-Stacking. Journal of the American Chemical Society 2002, 124, 12200–12209. (29) Spillebout, F.; B´egu´e, D.; Baraille, I.; Shaw, J. M. On Discerning Intermolecular and Intramolecular Vibrations in Experimental Acene Spectra. Energy & Fuels 2014, 28, 2933–2947. (30) Bree, A.; Zwarich, R. The Vibrations of Dibenzothiophene. Spectrochimica Acta Part A: Molecular Spectroscopy 1971, 27, 599–620. (31) Michaelian, K. H.; Billinghurst, B. E.; Shaw, J. M.; Lastovka, V. Far-Infrared Photoacoustic Spectra of Tetracene, Pentacene, Perylene and Pyrene. Vibrational Spectroscopy 2009, 49, 28–31. (32) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H. et al. Gaussian 09, Revision A.02. Gaussian, Inc., Wallingford CT 2016, 42

ACS Paragon Plus Environment

Page 42 of 48

Page 43 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(33) Chai, J.; Head-Gordon, M. Systematic Optimization of Long-Range Corrected Hybrid Density Functionals. The Journal of chemical physics 2008, 128, 084106. (34) Salzner, U.; Aydin, A. Improved Prediction of Properties of ⇡-Conjugated Oligomers with Range-Separated Hybrid Density Functionals. Journal of chemical theory and computation 2011, 7, 2568–2583. (35) Wiberg, K. B. Basis Set E↵ects on Calculated Geometries: 6-311++ G** vs. aug-ccpVDZ. Journal of computational chemistry 2004, 25, 1342–1346. (36) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of Van der Waals Complexes. Chemical Reviews 1994, 94, 1887–1930. (37) Hesselmann, A.; Jansen, G.; Sch¨ utz, M. Density-Functional Theory-SymmetryAdapted Intermolecular Perturbation Theory with Density Fitting: A New Efficient Method to Study Intermolecular Interaction Energies. The Journal of chemical physics 2005, 122, 014103. (38) Podeszwa, R.; Bukowski, R.; Szalewicz, K. Potential Energy Surface for the Benzene Dimer and Perturbational Analysis of ⇡-⇡ Interactions. The Journal of Physical Chemistry A 2006, 110, 10345–10354. (39) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Sch¨ utz, M. Molpro: A GeneralPurpose Quantum Chemistry Program Package. Wiley Interdisciplinary Reviews: Computational Molecular Science 2012, 2, 242–253. (40) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0 Model. The Journal of chemical physics 1999, 110, 6158–6170.

43

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(41) Barone, V.; Biczysko, M.; Bloino, J. Fully Anharmonic IR and Raman Spectra of Medium-Size Molecular Systems: Accuracy and Interpretation. Physical Chemistry Chemical Physics 2014, 16, 1759–1787. (42) Michaelian, K.; Oladepo, S.; Shaw, J.; Liu, X.; B´egu´e, D.; Baraille, I. Raman and Photoacoustic Infrared Spectra of Fluorene Derivatives: Experiment and Calculations. Vibrational Spectroscopy 2014, 74, 33–46. (43) Choi, S. E.; Light, J. C. Highly Excited Vibrational Eigenstates of Nonlinear Triatomic Molecules. Application to H2 O. The Journal of chemical physics 1992, 97, 7031–7054. (44) Bowman, J. M.; Gazdy, B. A Truncation/Recoupling Method for Basis Set Calculations of Eigenvalues and Eigenvectors. The Journal of chemical physics 1991, 94, 454–460. (45) Rauhut, G. Configuration Selection as a Route Towards Efficient Vibrational Configuration Interaction Calculations. The Journal of chemical physics 2007, 127, 184109. (46) Mizukami, W.; Tew, D. P. A Second-Order Multi-Reference Perturbation Method for Molecular Vibrations. The Journal of chemical physics 2013, 139, 194108. (47) Sibaev, M.; Crittenden, D. L. Balancing Accuracy and Efficiency in Selecting Vibrational Configuration Interaction Basis States Using Vibrational Perturbation Theory. The Journal of Chemical Physics 2016, 145, 064106. (48) Sibaev, M.; Crittenden, D. L. PyVCI: A Flexible Open-Source Code for Calculating Accurate Molecular Infrared Spectra. Computer Physics Communications 2016, 203, 290–297. (49) Lanczos, C. An Iteration Method for the Solution of the Eigenvalue Problem of Linear Di↵erential and Integral Operators; United States Governm. Press Office Los Angeles, CA, 1950.

44

ACS Paragon Plus Environment

Page 44 of 48

Page 45 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(50) Wyatt, R. E. Matrix Spectroscopy: Computation of Interior Eigenstates of Large Matrices Using Layered Iteration. Physical Review E 1995, 51, 3643. (51) Wang, X.-G.; Carrington Jr, T. Vibrational Energy Levels of CH5+. The Journal of chemical physics 2008, 129, 234102. (52) Seidler, P.; Hansen, M. B.; Gy˝or↵y, W.; To↵oli, D.; Christiansen, O. Vibrational Absorption Spectra Calculated fromVvibrational Configuration Interaction Response Theory Using the Lanczos Method. The Journal of chemical physics 2010, 132, 164105. (53) Avila, G.; Carrington Jr, T. Using Nonproduct Quadrature Grids to Solve the Vibrational Schr¨odinger Equation in 12D. The Journal of chemical physics 2011, 134, 054126. (54) Garnier, R.; Odunlami, M.; Le Bris, V.; B´egu´e, D.; Baraille, I.; Coulaud, O. Adaptive Vibrational Configuration Interaction (A-VCI): A Posteriori Error Estimation to Efficiently Compute Anharmonic IR Spectra. The Journal of chemical physics 2016, 144, 204123. (55) Odunlami, M.; Le Bris, V.; B´egu´e, D.; Baraille, I.; Coulaud, O. A-VCI: A Flexible Method to Efficiently Compute Vibrational Spectra. The Journal of Chemical Physics 2017, 146, 214108. (56) Hansen, M. B.; Sparta, M.; Seidler, P.; To↵oli, D.; Christiansen, O. New Formulation and Implementation of Vibrational self-Consistent Field Theory. Journal of chemical theory and computation 2009, 6, 235–248. (57) Leclerc, A.; Carrington, T. Calculating Vibrational Spectra with Sum of Product Basis Functions without Storing Full-Dimensional Vectors or Matrices. The Journal of chemical physics 2014, 140, 174111.

45

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(58) Barone, V. Anharmonic Vibrational Properties by a Fully Automated Second-Order Perturbative Approach. The Journal of chemical physics 2005, 122, 014108. (59) Begue, D.; Elissalde, S.; Pere, E.; Iratcabal, P.; Pouchan, C. New Theoretical and Experimental Infrared Results on Formaldehyde in Solution. The Journal of Physical Chemistry A 2006, 110, 7793–7800. (60) B´egu´e, D.; Gohaud, N.; Pouchan, C.; Cassam-Chena¨ı, P.; Li´evin, J. A Comparison of two Methods for Selecting Vibrational Configuration Interaction Spaces on a Heptatomic System: Ethylene Oxide. The Journal of chemical physics 2007, 127, 164115. (61) Wang, Y.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Full-Dimensional, ab initio Potential Energy and Dipole Moment Surfaces for Water. The Journal of chemical physics 2009, 131, 054511. (62) Baraille, I.; Larrieu, C.; Dargelos, A.; Chaillet, M. Calculation of Non-Fundamental IR Frequencies and Intensities at the Anharmonic Level. I. The Overtone, Combination and Di↵erence Bands of Diazomethane, H2 CN2 . Chemical Physics 2001, 273, 91–101. (63) Scribano, Y.; Benoit, D. M. Iterative Active-Space Selection for Vibrational Configuration Interaction Calculations Using a Reduced-Coupling VSCF Basis. Chemical Physics Letters 2008, 458, 384–387. (64) Brown, J.; Carrington Jr, T. Using an Expanding Nondirect Product Harmonic Basis with an Iterative Eigensolver to Compute Vibrational Energy Levels with as Many as Seven Atoms. The Journal of Chemical Physics 2016, 145, 144104. (65) Klots, T.; Collier, W. Heteroatom Derivatives of Indene. Part 2. Vibrational Spectra of Benzothiophene and Benzothiazole. Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 1995, 51, 1273–1290.

46

ACS Paragon Plus Environment

Page 46 of 48

Page 47 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(66) Sodero, A. C.; Santos Silva, H.; Guevara Level, P.; Bouyssiere, B.; Korb, J.-P.; Carrier, H.; Alfarra, A.; Begue, D.; Baraille, I. Investigation of the E↵ect of Sulfur Heteroatom on Asphaltene Aggregation. Energy & Fuels 2016, 30, 4758–4766. (67) Silva, H. S.; Sodero, A. C.; Bouyssiere, B.; Carrier, H.; Korb, J.; Alfarra, A.; Vallverdu, G.; B´egu´e, D.; Baraille, I. Molecular Dynamics Study of Nanoaggregation in Asphaltene Mixtures: E↵ects of the N, O, and S Heteroatoms. Energy & Fuels 2016, 30, 5656–5664. (68) Johnson, E. R.; Keinan, S.; Mori-Sanchez, P.; Contreras-Garcia, J.; Cohen, A. J.; Yang, W. Revealing Noncovalent Interactions. Journal of the American Chemical Society 2010, 132, 6498–6506.

47

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry

Structural parameters of benzothiophene dimers.

48

ACS Paragon Plus Environment

Page 48 of 48