Article pubs.acs.org/JPCA
Possible Formation of Metastable PAH Dimers upon Pickup by Helium Droplets F. Calvo,*,†,‡ E. Yurtsever,§ and Ö . Birer§ †
Univ. Grenoble Alpes, LIPHY, F-38000 Grenoble, France CNRS, LIPHY, F-38000 Grenoble, France § Koç University, Chemistry Department, Rumeli Feneri Yolu, Sariyer 34450, Istanbul, Turkey ‡
S Supporting Information *
ABSTRACT: Using path-integral molecular dynamics simulations and two quantum-mechanical-based force fields, we have investigated the conformational stability of dimers of a polycyclic aromatic hydrocarbon, perylene (C20H12), produced under typical experimental conditions of successive pick-up under helium nanodroplet environment. The most stable configurations are found to be of the stacked form with different relative orientations of the main molecular axes, perpendicular or T-shaped dimers being energetically much disfavored; however, in the presence of helium our simulations suggest that the time for rearrangement and π-stacking may be rather long and exceed hundreds of picoseconds. In addition, highly metastable dimers that are stacked but with a helium monolayer sandwiched between the two molecules are also found as likely products upon successive pickup. This stabilization occurs owing to the stronger localization of the helium atoms facing the aromatic rings, which is further enhanced in the dimer. The implications of the present results are discussed in the perspective of possible identification by spectroscopic methods.
1. INTRODUCTION Helium droplets are an ideal, chemically inert solvent to study molecular systems near their ground rovibrational state. Because the technique was introduced in the early 1990s,1 it has provided a wealth of information into the high-resolution spectroscopy of solute molecules or ions ranging from atomic impurities1−6 to clusters of fullerenes.7 The phenomenon of superfluidity, in particular, has been scrutinized in great detail through its spectroscopic influence8,9 and its possible partial suppression in the vicinity of ionic dopants, leading to the formation of so-called snowballs.10−12 Several molecular complexes have been previously investigated under helium nanodroplet environment. Besides the aforementioned rather extreme case of fullerene clusters, hydrogen cyanide13 was found to form linear clusters bound by dipole−dipole interactions, whereas the onset of acidic behavior in stepwise hydrated HCl was studied both experimentally and theoretically14−17 and shown to occur quite early with only four water molecules. Such complexes are formed by successive pick-up of individual molecules, with the initially pure helium droplet passing through a gas chamber containing the solute molecules. The excess energy of the capture process is released by the evaporation of solvent atoms, with the temperature being thus maintained below a few Kelvins. The complex is subsequently formed by the aggregation of the molecules together that are bound by long-range forces. © XXXX American Chemical Society
In small complexes, the potential energy surfaces are usually rather simple and the molecular complexes are expected to lie in their lowest energy structure; however, as the number of molecules increases the energy landscape becomes increasingly rugged,18 and the complexes are more likely to be trapped into higher energy conformations, with the time needed for relaxation to the global minimum being excessively long under the cryogenic temperatures of the helium solvent. The formation of energetically metastable conformations was already observed by Nauta and Miller in the case of hydrogen cyanide13 and interpreted as being caused by the long-range interactions between these highly dipolar molecules that drive the aggregation toward the entropically favored linear oligomers. Recently, an experiment was reported in which the perylene dimer could be formed in helium droplets and spectroscopically characterized,19 with the results being compared with ab initio calculations in the gas phase. Support was found for the predominance of parallel displaced stacked geometries, although the contribution of the higher-lying rotated stacked conformer could not be ruled out with the available resolution of both calculation and experiment.19 Those results contrast with previous conclusions reached for the related peryleneReceived: December 18, 2015 Revised: February 16, 2016
A
DOI: 10.1021/acs.jpca.5b12394 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
Here we consider the perylene dimer as a set of (twice) N = 32 carbon and hydrogen atoms with classical configuration collectively denoted as RHC,i, i = 1 or 2, surrounded by a set of nHe = 500 4He atoms located at position vector rHe = {rHe,i}, with rHe,i being the classical position of helium atom i. For brevity, we further denote by R = {RHC,i, rHe,j} the entire set of positions of all 2N + nHe atoms in the system. 2.1. Potential Energy Surface. The total interaction energy of the system at classical collective position R is written in the additive approximation as
3,4,9,10-tetracarboxylic dianhydride (PTDCA) complexes, which were suggested to adopt T-shape conformers owing to dipolar side groups.20 The present work aims to theoretically evaluate the influence of the helium solvent on the possible conformations adopted by such a molecular complex, taking into account the pick-up nature of the formation process, that is, not assuming from the start that the complex readily finds its global energy minimum. In a realistic pick-up experiment, the different monomers are captured at the surface of the droplets and have to make their way to aggregate to each other at the center. While the interactions with helium themselves are expected to be small enough not to alter the energy landscape of the complex, we show that this view is misleading because the solvent could still contribute to forming dimers by slowly relaxing transient structures on the way to the stable π-stacked minima. Furthermore, our simulations indicate that metastable stacked conformers with helium atoms sandwiched between two monomers may be produced owing to the stronger vibrational localization of helium in the vicinity of aromatic rings. The strengthening of helium localization near aromatic rings was noted decades ago with the reentrant freezing of helium monolayers on graphite,21 and, in the context of helium droplets, on the stronger localization and suppression of superfluidity of 4He around various aromatic molecules such as benzene,22,23 anthracene,24 phthalocyanine,25 and much larger polycyclic aromatic hydrocarbons (PAHs)26,27 ranging up to circumcoronene28 as well as around 3D fullerene molecules.29,30 The increased localization of helium atoms on PAH molecules has been evidenced in hole-burning experiments through line splittings.31,32 The approach employed to address those issues is based on the framework of path-integral molecular dynamics (PIMD) and more particularly ring-polymer MD, which captures the actual time-dependent processes to some extent.33,34 Although our objectives are not to quantify in detail the statistical occurences of the various conformers or their formation rates in the pick-up processes, the PIMD approach is practical because it can deal with realistic droplet sizes of several hundreds of atoms, over time scales exceeding tens of picoseconds. One important limitation, however, is that we ignore exchange statistics here, as this would require much more sophisticated techniques35,36 involving permutation sampling that have the drawback of making kinetic processes harder to address. The article is organized as follows. In the next section, we present the computational methods used to simulate the perylene dimer under helium environment in the normal fluid state, including the fully atomistic potential energy surface and some algorithmic details of the PIMD approach. The tools necessary for geometric and spectroscopic characterizations are also laid out in this section. The main results are presented and discussed in Section 3, focusing on different conformers formed through diverse mechanisms likely to all be present in pick-up processes. Finally, we summarize and conclude in Section 3 about the possible relevance of the present results to other types of molecular complexes and their possible experimental detection through spectroscopic methods.
V (R) =
∑ VTB(R HC,i) + i
+
∑
Vinter(rij)
i ∈ HC,1; j ∈ HC,2
∑
VHe(rij) +
i , j ∈ He, i < j
∑ i ∈ HC, j ∈ He
Vpair(rij) (1)
where the successive terms account for the intramolecular and intermolecular interactions within the hydrocarbons, for the helium−helium and hydrocarbon−helium interactions, respectively. The intramolecular interactions VTB within each perylene molecule are taken here using the tight-binding (TB) quantum model developed by Van-Oanh et al.,37 which provides a rather accurate description of the spectroscopic properties of polyaromatics. For the intermolecular pairwise contributions Vinter, and following previous work on PAH clusters,38 we assume that each atom carries a Lennard-Jones center and a partial charge representing the multipolar electrostatic distribution. The LJ contribution accounts for the long-range dispersion attraction and the short-range Pauli repulsion, while the multipolar distribution plays a role in stabilizing specific orientations between the two molecules near equilibrium. Following a standard procedure, these partial charges were obtained from a DFT calculation at the B3LYP/6-31+G* level using the Gaussian09 software package39 and the RESP definition. They are provided as Supporting Information together with the monomer geometry at equilibrium. The Lennard-Jones intermolecular forces between perylene molecules are not straightforward to set or even parametrize due to the difficulty of evaluating the long-range interactions in such large systems and partition them into atomic contributions. Here, as a way to further provide some rough estimator for the robustness of our findings, all simulations were performed with the same set of LJ parameters as in ref 38, as previously used by van de Waals for benzene clusters,40 but were repeated using the OPLS set of LJ parameters41 otherwise under the same conditions. The interactions between helium and hydrocarbon atoms were taken following our recent work on helium-coated PAH cations28 without further modification, and the Janzen−Aziz potential42 was used to model the interaction among helium atoms. Because of their very weak interactions, helium atoms tend to easily evaporate even in small droplets. This causes convergence issues in the thermostated simulations and was fixed by imposing an additional but soft confining potential Vrep to the global interaction energy, acting only beyond a sufficiently large radius to not exert pressure on the system of interest. For each atom, Vrep acts on the radial distance ri of helium atom i to the origin of the reference frame accordingly with
2. METHODS Our approach to the formation of molecular dimers in helium environment relies on the path-integral molecular dynamics framework using explicit potential energy surfaces with atomistic force fields.
Vrep(ri) = κ(ri − R 0)4 B
(2) DOI: 10.1021/acs.jpca.5b12394 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A where R0 = 20 Å the radius of the container and κ = 10−5 atomic units. 2.2. Path-Integral Molecular Dynamics. Path-integral molecular dynamics is a rigorous method that incorporates quantum effects such as delocalization and tunneling in the description of a many-body system in canonical equilibrium. After the technique was introduced in molecular simulation,43,44 several related approaches were subsequently proposed, aiming to evaluate time-dependent properties from quantum time autocorrelation functions, most notably the centroid MD (CMD) method of Voth and coworkers45 and the ring-polymer MD (RPMD) method of Manolopoulos and coworkers.33 It is not our purpose here to review those methods46 but only to give the main justification of using PIMD in the first place. While path-integral Monte Carlo is probably the only practical approach to study those systems at thermal equilibrium, it lacks any dynamical component and would be less convenient to address possibly metastable configurations. Although we are aware that regular (canonical) PIMD does not strictly provide any information about the rates at which molecular processes take place, its CMD and RPMD versions are designed to approach the true dynamics.33,34,45,47 Here we have chosen the RPMD method, which we have successfully used in the past for a variety of other problems involving polyaromatics29,48 to get insight into the time evolution of molecular configurations initially chosen from different families typically occurring along successive pick-up processes. The limitations of both CMD and RPMD for producing quantitative results are well-documented,49 and our ambitions are not to calculate rates but at best to evaluate approximate time scales for possible rearrangements. In the path-integral framework, the n-atom molecular system is replaced by a higher dimensional n × P set of particles, with the various replicas interacting successively through harmonic bonds on each atom and through the physical forces within each replica. In brief, both equilibrated PIMD and ring-polymer MD methods amount to solving the classical equations of motion for this extended dimensional system, either in the presence of (massive) thermostats in PIMD or having selected initial conditions that are already thermally equilibrated in RPMD. Here we further take benefit of the additive nature of the interaction energy of eq 1 by employing the ring contraction technique50 to describe the different degrees of delocalization of the perylene dimer and the helium environment, choosing P = 32 for the former and P = 64 for the latter at the temperature T = 1 K chosen in this study. While such values of P are probably not large enough to quantify energetic properties, considering the rather large size of the system under study, they are sufficient to describe the structural properties we are mainly interested in here. The simulation protocol was as follows for all trajectories. The two perylene molecules were initially placed relative to each other either in the most stable stacked structure or on their way to equilibration, assuming different conformations such as T shapes or parallel stacks but separated by a larger distance exceeding 6 Å. The perylene molecules were then placed into a helium bath, borrowing the helium centroid positions from an external independent PIMD trajectory on pure 500-atom 4He clusters. Helium atoms too close to the perylene molecules were removed and placed at the outer regions within the container border, R0. A short thermostated
PIMD trajectory was then initiated and propagated for 1 ps before the thermostat was released and the RPMD trajectory was further propagated. Because of the significant computational cost of integrating the PIMD and RPMD equations of motion for such large-size systems, only limited numbers of trajectories could be performed, although generally producing results that are consistent among them. 2.3. Electronic Spectra. All quantum calculations were carried out with Gaussian09 package.39 The optical spectra of monomer, dimers and dimers in the presence of helium atoms were obtained from TD-DFT using the M06-2X exchange functional and the cc-pVDZ basis set. Noncovalent interactions are usualy problematic within DFT methodology; however, M06-2X was shown to be reliable in describing such long-range interactions.51 Furthermore, for a set of benchmark molecules, it was shown that the most accurate optical spectra were obtained from TD-DFT and M06-2X functional combination.53 2.4. Computational Details. All path-integral and RPMD trajectories employed a time step of 0.2 fs and were propagated over 50 ps each. They were mostly analyzed in terms of geometries, focusing on the possible rearrangements of the perylene dimer rather than on the response of helium to specific initial conformations, even though the helium atoms turn out to participate actively in the stabilization of some of these conformations (vide infra). In particular, at the very low temperatures imposed by the helium environment, the individual hydrocarbon molecules are essentially frozen and vibrate only through zero-point contributions; however, the intermolecular interactions are much softer, and it is reasonable to describe the dimer conformations from the centroids only, treating them as classical in the structural analysis. In the simulations, we associate an orthonormal reference frame to each molecule from the centroid positions, choosing the unit vectors ux⃗ and vx⃗ of the main molecular axes as joining the centers of mass (also center of frame) to specific outer carbon atoms, as indicated in Figure 1. The vectors uy⃗ and vy⃗ defining the short axes and the molecular planes are obtained from the midbond carbon dimer roughly parallel to the long axis vectors, and finally the vectors
Figure 1. Reference frame used to identify the relative motion of the molecules in the perylene dimer: axes vectors associated with the two molecules and vector joining their centers of mass. C
DOI: 10.1021/acs.jpca.5b12394 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A u⃗z and vz⃗ complete the orthonormal frame. To monitor the relative motion between the two monomers, we also required the vector r1⃗ 2 joining the two centers of mass. From the two orthonormal frames several angles between unit vectors are defined as θx = cos−1 u⃗x·vx⃗ , with θy and θz being obtained similarly. The shearing motion of the two monomers under parallel stacked conformation can be evidenced from the angle between either of the normal vectors u⃗z or vz⃗ with the (normalized) intermolecular vector r̃12 = r 1⃗ 2 /||r 1⃗ 2 || as θr = cos−1 u⃗z·r̃12. These geometric parameters were evaluated as a function of time for all RPMD trajectories based on the centroids of the perylene molecules. Additional quantities such as the helium density were also determined for specific conformations for which a common reference frame can be conveniently defined.
other, for (perpendicular) T-shaped dimers, or for intermediate “Y-shaped” conformations. Changing the intermolecular interactions between the two LJ sets does not drastically alter the energy landscapes of the perylene dimer; however, the relative stabilities of the two conformers are somewhat affected, as shown in Table 1, where their energies and geometrical properties are given in terms of the angles θx, θy, and θr. While the two isomers are nearly degenerate with the vdW model, the parallel displaced conformation lies 55 cm−1 above the rotated stack with the OPLS set of LJ parameters, with the electrostatic interactions remaining unchanged. In both cases we expect the energy surface to be very flat parallel to the molecules, and rough estimates of the energy barrier connecting them based on linear interpolation leads to values below 100 cm−1 for both LJ models. The geometric properties of the two isomers barely depend on the intermolecular forces, and as can be seen from Table 1, the three angles θx, θz, and θr sufficiently differ for the conformers to be distinguishable from one another. The very flat nature of the energy surface can be gauged by considering the dynamics of the perylene dimer in vacuum. Initiating each PIMD trajectory at T = 1 K from the corresponding most stable structure, the variations of the three angles are depicted as a function of time in Figure 3 for the two models. Each angle exhibits fluctuations that are quite significant over time, except θz, which remains lower than 5° in general. The low magnitude of this angle shows that the molecules remain stacked during their relative motion, which is essentially lateral. The much larger fluctuations obtained for both θx and θr span the two values characteristic of the parallel
3. RESULTS AND DISCUSSION The present potentials predict only two stable conformations for the perylene dimer in vacuum. The binding energies relative to two monomers at infinite separation are 0.74 and 0.56 eV with the vdW and OPLS intermolecular models, respectively. (See Table 1.) Both dimers are stacked and of the parallel Table 1. Relative Energies and Geometric Properties of the Two Locally Stable Stacked Conformations of the Perylene Dimer Found with the Present Potential Energy Surfacesa property
parallel displaced
rotated stack
energy (vdW model) energy (OPLS model) θx θz θr
+2.7 cm−1 0.56 eV 57.3° 1.9° 15.1°
0.74 eV +55 cm−1 40.1° 0.2° 0.2°
a
Energies in electronvolts denote the binding energies in the most stable conformation relative to infinite separation, and for each model the energy in cm−1 is that of the higher-lying conformation relative to the most stable one. The angles show minor variations with the LJ model.
displaced or rotated types previously identified using MP2 electronic structure calculations.19 They are displayed in Figure 2. No evidence was found from classical or PIMD simulations for stable stacked dimers with aromatic rings on top of each
Figure 3. Time variations of the angles θx, θy, and θr (see text for their definitions) as obtained from RPMD trajectories of the isolated perylene dimer at T = 1 K using the two intermolecular LJ force fields with parameters vdW or OPLS. The corresponding values of those angles in the stable structure are denoted by PD and R on the right for parallel displaced and rotated stacked structures, respectively.
Figure 2. Two stable stacked forms of the perylene dimer in vacuum predicted by the models. The long and short axes of each molecule are displayed to highlight the relative locations of the centers of mass and the orientations. D
DOI: 10.1021/acs.jpca.5b12394 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
monitored from the angle θz alone, which is sufficient to measure the extent of stacking. Four independent RPMD trajectories were performed for each LJ model, the results of which are displayed in Figure 5 with the angle θz as a function
displaced and rotated stacked conformers, indicating that both are visited across the trajectories and for both models. Visual inspection confirms this analysis as well as systematic quenching of the centroid positions (results not shown). Placing now the stacked dimers into a 500-atom droplet of 4 He at the same temperature of 1 K and repeating the PIMD simulations, the geometrical parameters exhibit much lower fluctuations, as depicted in Figure 4. In the presence of helium
Figure 5. Time rearrangement of a perylene dimer initially having a Tshape configuration (c) surrounded by 500 He atoms, as monitored from the θz angle along four representative PIMD trajectories with the two intermolecular LJ force fields with parameters vdW (a) and OPLS (b). The final configurations are of the types depicted in (d) bent shape or (e) stacked form.
Figure 4. Same as Figure 3 but for stacked perylene dimers surrounded by 500 4He atoms.
of time. Under such initial conditions, the two intermolecular models predict rather different kinetics, although the long-time asymptotic behavior should be similar, with the eventual stacking of the dimer. With the vdW model the molecules rotate within
