MM Approach to Isomerization of Ruthenium(II) Sulfur Dioxide

Subscriber access provided by UNIV OF ALABAMA BIRMINGHAM

C: Surfaces, Interfaces, Porous Materials, and Catalysis

QM/MM Approach to Isomerization of Ruthenium(II) Sulfur Dioxide Complex in Crystal; Comparison with Solution and Gas Phases Shinji Aono, and Shigeyoshi Sakaki J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b04774 • Publication Date (Web): 14 Aug 2018 Downloaded from http://pubs.acs.org on August 19, 2018

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 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 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.

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 55 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

QM/MM Approach to Isomerization of Ruthenium(II) Sulfur Dioxide Complex in Crystal; Comparison with Solution and Gas Phases Shinji Aono and Shigeyoshi Sakaki∗ Fukui Institute for Fundamental Chemistry, Kyoto University, Nishihiraki-cho, Takano, Sakyo-ku, Kyoto 606-8103, Japan E-mail: [email protected]

∗ To

whom correspondence should be addressed

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 Theoretical study of chemical reaction in molecular crystal is challenging research target. To optimized transition state in crystal, we improved crystal model and modified quantum mechanics/molecular mechanics (QM/MM) method based on periodic MM crystal model. We applied this method to ruthenium(II) sulfur dioxide complexRuthenium(II) sulfur dioxide complex [RuII (NH3 )4 (SO2 )Py]2+ 1Py which is potentially useful for optical data storage system in crystal because this Ru complex has two isolable meta-stable structures, η 1 -O-bound MS1 and η 2 -S,O-bound MS2, in addition to η 1 -S-bound ground state GS in crystal. MS2 and MS1 are formed from GS by photo-irradiation at low temperature (10 - 100 K) and thermal isomerization occurs from MS2 to MS1 and from MS2 to GS by raising temperature. For comparison, the thermal isomerization of [RuII (NH3 )4 (SO2 )L]n+ (L = Cl− , W, and Py) was investigated in gas and aqueous phases, where three-regions 3D reference site interaction model self-consistent field (RISM-SCF) method was employed to incorporate solvation effects in aqueous phase. In gas phase, MS2 is more stable than MS1 for L = Cl− but less stable for L = W and Py. Thermal isomerization occurs with similar changes in geometry and charge distribution among L = Cl− , W, and Py, whereas energy profile of L = Cl− differs very much from those of L = W and Py. In aqueous phase, MS1 becomes much less stable than MS2 and the RuII complex becomes lower in energy following the order MS1 > MS2 > GS for L = Cl− , W, and Py because of small solvation free energy for MS1. In crystal, 1Py also becomes lower in energy following the order MS1 > MS2 > GS. The MS1→MS2 isomerization occurs with much smaller activation energy than the MS2→GS in crystal phase, which is consistent with the experimental observation. In crystal, short-range steric repulsion between SO2 of target QM molecule and neighboring Ru complex and long-range ES interaction between the target molecule and MM crystal play crucially important roles in determining relative stabilities of MS1 and MS2 and energy barrier of the thermal isomerization.

2


Page 2 of 55

Page 3 of 55 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


Introduction Recently, spin state, absorption and emission spectra in solid system such as crystal and amorphous systems receive a lot of attentions in both of experimental and theoretical chemistry because these properties are utilized for molecular devices and functional materials. In many systems, molecular property and reactivity in solid system significantly differ from gas and solution phases, because molecule in solid system is much more strongly influenced by surrounding molecules through long-range electrostatic (ES) and short-range steric interactions than in gas and solution phases. In this regard, correct evaluation of surrounding molecular effects is necessary for understanding, utilizing, and designing molecules in solid system. Photo-induced linkage isomerization, in which a ligand bound to a metal center has one or more meta-stable structures different from equilibrium one, has been experimentally and theoretically investigated because such isomerization is potentially useful for optical data storage, 1–3 if ground state of isomer can be assigned as “0" and the meta-stable one as “1". In recent years, metastable structures of such materials have been experimentally determined by development of photocrystallographic technique, 4–6 in which various active ligands such as NO− , 5,7–13 NO2 , 7,14–16 N2 , 17,18 dimethyl sulfoxide, 19,20 and other organic sulfoxides 21–23 behave as photo-active linkage isomer. Ruthenium(II) sulfur dioxide complex 1, [RuII (NH3 )4 (SO2 )L]n+ (n = 1 for L = Cl− ; n = 2 for L = H2 O (W) and pyridine (Py)), is one of the most interesting systems exhibiting photo-induced linkage isomerization in crystal, because the SO2 moiety has several coordination structures, as follows: 24–30 Ru(II) ion with d6 electron configuration has unoccupied dz2 and dxy orbitals which are useful to form coordination bonds with SO2 and L on the z-axis and four NH3 ligands on the xy-plane. It has η 1 -S-bound SO2 coordination structure at ground state (GS), as shown in Scheme 1. In addition, these Ru complexes have two isolable photo-induced meta-stable structures such as

η 1 -O-bound MS1 structure and η 2 -S,O-bound MS2 structure. These isomers are possible because the SO2 moiety has several frontier orbitals; the HOMO (lσ ) is a lone pair type orbital containing

σ ∗ -antibonding overlap between S and O atoms, which participates in σ -donation from S of SO2 3



to Ru dz2 in GS. The next HOMO is in-plain non-bonding p orbital (nin ), which participates in ∗ orbital, σ -donation from O of SO2 to Ru dz2 in MS1. The LUMO is out-of-plane antibonding πout

which participates in π -back donation from doubly occupied Ru dπ orbital to SO2 π ∗ orbital in MS2, as shown in Scheme 2. At low temperature (10 - 100 K), GS is converted to MS2 and MS1 by photo-irradiation and then thermal isomerization occurs from MS1 to MS2 and from MS2 to GS by raising temperature; see Scheme 1. [RuII (NH3 )4 (SO2 )(Py)]2+ (2Cl− ·H2 O) in crystal has been experimentally reported to exhibit large quantum yield of MS1 by visible light irradiation (425 - 650 nm) at 100 K; this complex is considered the best candidate for application because much lower temperature than 100 K is needed to observe experimentally MS1 in other RuII complexes. 29 GS, MS1, and MS2 were observed to be 45, 40, and 15 %, respectively, after photo-conversion. When temperature was elevated from 100 to 140 K, the amount of GS changed little but MS1 disappeared and that of MS2 largely increased to 50 % after several hours, indicating that MS1 decays to MS2 in 100 140 K but MS2 does not to GS in this temperature. When temperature was elevated from 150 to 200 K, MS2 decayed to GS. The isomerization behavior is largely influenced by L and Y of [RuII (NH3 )4 (SO2 )L]Y and crystal structure. Hence, the surrounding molecular effects must be taken into consideration to understand this photo-induced linkage isomerization in crystal. However, theoretical study of such chemical reaction in crystal has not been performed yet by incorporating sufficiently surrounding molecular effects. 24,27,28 In this study, the above-mentioned thermal isomerization of [RuII (NH3 )4 (SO2 )Py]2+ (2Cl− · H2 O) in crystal was theoretically investigated, using quantum mechanics/molecular mechanics (QM/MM) method based on periodic MM crystal model with self-consistent point charge distribution. 31,32 In this method, crystalline effects can be effectively incorporated in calculation by embedding QM molecule in periodic MM crystal, where point charges and geometries of MM crystal are determined in the self-consistent filed (SCF) manner by applying the translational and symmetrizing operations to QM molecule in minimum unit cell. We also investigated thermal isomerization of [RuII (NH3 )4 (SO2 )L]n+ (L = Cl− , W, and Py) in gas and aqueous phases to elu-

4


Page 4 of 55

Page 5 of 55 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


cidate how much surrounding molecules influence thermal isomerization in crystal and solution phases and disclose characteristic feature of reaction in gas, aqueous, and crystal phases, where three-regions 3D-RISM-SCF method 33–35 was employed to incorporate solvation effects in aqueous phase.

Theoretical Method Since the basic concept of QM/MM method was introduced in 1976 by Warshel and Levitt, 31 the QM/MM approach has been established as a powerful tool for modeling large biomolecular systems, chemical reactions, and molecular properties of large inorganic and organometallic compounds in explicit solvent system. The most important advantage of the QM/MM approach is its efficiency to study the kinetic and entropic effects of proteins and solution, because MM simulation can be performed efficiently (scale of O(N) to O(N 2 )) with accurate description of QM region. In the present study, we applied the QM/MM method to chemical reaction in crystal. The purpose of the QM/MM method is to incorporate static effects of long-range electrostatic interaction with crystal and short-range steric effects by neighboring molecules on the target QM molecule, which is a bit different from the usual QM/MM method applied to biomolecular and solution systems.

Construction of MM crystal model To incorporate crystalline effects on reaction of QM molecule, we constructed MM crystal model in closed-shell singlet ground state by using the periodic QM/MM calculation method recently proposed by us, 32 where self-consistent point charge was employed to describe MM region. Though the same idea was independently proposed by Iwano and Shimoi, 36–38 we reported the algorithm to perform the geometry optimization with energy gradient considering the derivatives of selfconsistent point charge. We will briefly explain here the minimum of this method, because some modification of crystal model was needed here to optimize the transition state of chemical reaction in crystal.

5



Page 6 of 55

The key idea is based on the assumption that all the “symmetrically minimal subunits" have the same charge distribution and geometry in molecular crystal, which is used to construct the minimum unit cell by applying the symmetrizing operation T sym . In other words, the point charges and coordinates {qa , ra } of all the molecules in crystal can be produced by those {qa′ , ra′ } of a given minimal subunit, using proper symmetrizing matrix Tsym , center of symmetry x0 , and/or translational vector Ttrans , when atom a is symmetrically and periodically equivalent to atom a′ in molecular crystal, as shown in Eq. (1); {qa , ra } = {qa′ , Tsym (r′a′ − x0 ) + Ttrans }.

(1)

If the intermolecular interaction energy in crystal can be represented by the sum of two-body I interactions, the potential energy of total crystal is defined as the sum of QM energy EQM of IJ symmetrically minimal subunit I and a half of QM-MM interaction energy ∑ EInt(QM-MM) between

the minimal subunit I and all the remaining subunits J, as shown in Eq. (2); total I I Ecrystal ∝ Ecrystal ≡ EQM +

1 IJ . EInt(QM-MM) 2 J(̸∑ =I)

(2)

I IJ In this method, EQM and EInt(QM-MM) are represented by Eqs. (3) and (4);

I EQM = ⟨ΦI | Hˆ gas | ΦI ⟩, ) ( √ IJ EInt(QM-MM) = ∑ qa′VaI′ + εa′UaI′ .

(3) (4)

a′ ∈J(̸=I)

Here, Hˆ gas is Hamiltonian in gas phase, ΦI is wave function of symmetrically minimal subunit I in crystal, VaI′ is electrostatic potential (ESP) from subunit I acting on site a′ , UaI′ is Lennard-Jones (LJ) potential from the subunit I acting on the site a′ , εa′ is LJ parameter of site a′ , where the site a′ belongs to subunit J. The first term of Eq. (4) is directly evaluated by one-electron integral in inner region near the QM molecule and approximately evaluated as the sum of Coulombic interactions between point-charges using the Ewald summation method in outer region which is distant from 6


Page 7 of 55 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 QM molecule. In the intermediate region between these two regions, “VaI′ " was smoothly connected by using the switching function; w(r) = 1 − x3 {10 − x(15 − 6x)} for 0 ≤ x ≡

| r − rcen | −rinner ≤ 1, router − rinner

(5)

where rcen is center of the QM molecule and rinner and router are the radius of the inner and outer spherical regions, respectively. 32 The essentially the same method was used in 3D-RISM and 3DRISM-SCF methods proposed for solvation effects; 34 see pages S2-S3 in ESI for brief history of RISM-SCF method. I In Eq. (3), ΦI was determined by minimizing the energy Ecrystal with variational principle,

where qa′ depends on ΦI (Eq. (1)) and hence ΦI and qa′ are simultaneously solved in an SCF manner. The geometry of symmetrically minimal subunit I was further optimized by using analytical I energy gradient of Ecrystal . 32 This computational scheme is shown in Scheme 3a.

Geometry optimization of target QM molecule in the presence of MM crystal The transition state structure must be optimized for one QM molecule but MM crystal model was fixed to be the same as that of the reactant, because in geometry of only one molecule changes in this isomerization; if geometries of all molecules change at one time, phase transition occurs. In this regard, transition state is different between MS1→MS2 isomerization and MS2→MS1 isomerization: In the MS1→MS2 isomerization, the MM crystal structure is taken to be the same as that of MS1, because all molecules take MS1 structure. Hence, the transition state optimization was performed without SCF procedures between QM and MM geometries and between QM and MM charges; see pages S4-S5 in ESI for details of the application to the isomerization in crystal. For such investigation of transition state, we must make some rearrangement, as follows: After constructing MM crystal which is surrounding target molecule, we optimized transition state using QM/MM method by taking the target molecule as QM region under the conditions that the structure and charge distribution of MM crystal model are frozen, as shown in Scheme 3b. 40 In this step,

7



Page 8 of 55

the geometry and electronic structure differ between the target QM molecule and the remaining molecules in the MM crystal model and hence the potential energy of total crystal does not satisfy Eq. (2). 39 The potential energy of total crystal can be represented by Eq. (6) in the same manner as the usual QM/MM method;

total Ecrystal = EQM + EInt(QM-MM) + EMM ,

(6)

where EQM is the sum of the self energy ESelf(QM) of the QM molecule and the interaction energy EInt(QM-QM) between QM molecules, EInt(QM-MM) is the sum of interaction energies between QM and MM molecules, and EMM is the sum of self energy of MM molecule and interaction energy between MM molecules. Under the above-mentioned conditions, we can construct the Fock matrix and optimize the geometry by minimizing the sum of the EQM and EInt(QM-MM) energies in Eq. (6), because EMM term is constant (MM crystal is frozen). EMM is hereafter taken as zero value and the superscript total is omitted for simplicity. Using the variational principle, the Fock matrix F “total" of Ecrystal µν

can be represented by adding one-electron operator term of the QM-MM interaction to the Fock gas

matrix in gas phase Fµν , as shown in Eq. (7); { gas

Fµν = Fµν −

∑

qMM b

b∈MM

Qµν ,a′ w(rb )Aµν (rb ) − (1 − w(rb )) ∑ | rb − ra′ | a′ ∈QM

} ,

(7)

where the last term is evaluated with the Ewald summation, Qµν ,a′ is charge operator generating QM point-charge on the site a′ , w(r) denotes switching function to smoothly connect ESP between the inner and outer regions, 34,35 and Aµν is one-electron integral represented by atomic basis functions χµ and χν ; see Eq. (8); ∫

Aµν (r) =

dr′ χµ∗ (r′ )

8

1 χν (r′ ). | r − r′ |


(8)

Page 9 of 55 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 analytical gradient of Ecrystal can be derived by Eq. (9); ] [ ∂ Ecrystal √ ∂ UbLJ ′ ∂ EQM ∂ Z a = + ∑ + qMM εb ∑ | rb − ra′ | b ∂ ra ∂ ra ∂ r ∂ ra a′ ∈QM a b∈MM }] { [ ′ Q ∂ µν ,a − ∑ Dµν ∑ qMM ,(9) w(rb )Aµν (rb ) − (1 − w(rb )) ∑ b ∂ ra | rb − ra′ | µν b∈MM a′ ∈QM where the first term is evaluated in the same manner as energy gradient in gas phase, Dµν is oneelectron density matrix, and Zb is nuclear charge on the site b.

Models and Computational Details The first step is to construct periodic MM crystal models, in which all the molecules take the same GS, MS1, or MS2 structure. This condition is necessary to model the periodic MM crystal without determining potential energy functions of bond lengths, bond angles, and torsion angles of MM molecules. For this calculation, the target Ru complex [RuII (NH3 )4 (SO2 )Py]2+ 1Py , the nearest neighboring one 1Py on the +z-axis, and two pairs of counterions and crystal water (2Cl− ·H2 O) were taken as QM part to describe the interaction between the SO2 moiety of the target Ru complex and the pyridine moiety of neighboring Ru complex at the QM level, as shown in Scheme 4a. Then, we investigated the thermal isomerization of the target Ru complex 1Py using the MM crystal determined above. Transition state geometry of the target QM molecule was optimized, employing the target 1Py (QM0), the nearest neighboring two 1Py (QM1), and five Cl− counterions and one crystal H2 O (QM2) as QM part of the QM/MM calculation, as shown in Scheme 4b. This larger QM model was used to evaluate well the interaction of the target QM molecule with neighboring Ru complexes at the QM level. In this optimization, the periodic MM crystal constructed in the first step was frozen; 40 note that the relaxation of molecules surrounding QM0 is considered by optimizing the neighboring molecules QM1 and QM2. 41 The CCSD(T) potential energy of this trimer was evaluated by Eq. (10) using the ONIOM procedure, 42,43 where the energy of target QM molecules and the interaction energy between two QM molecules are calculated 9



Page 10 of 55

at the CCSD(T) and DFT levels, respectively. CCSD(T)

Ecrystal

B3PW91 (QM0 + QM1 + QM2) ≈ Ecrystal (QM0 + QM1 + QM2) + ( ) CCSD(T) B3PW91 Ecrystal (QM0) − Ecrystal (QM0) .

(10)

We implemented our own in-house code of the QM/MM method into GAMESS-US (version 1-MAY-2013) program. 44 Hessian matrix was updated by Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm for the geometry optimization of stable and meta-stable geometries. 45 POWELL’s quadratically convergent algorithm was employed for optimization of transition state. 45 Lattice vectors and initial structure of molecular crystal were taken from experimental data observed at 100, 120, and 140 K. 29,46 For QM calculation, we employed DFT method with the B3PW91 hybrid functional which consists of the Becke’s three-parameter exchange functionals 47 and the Perdew-Wang 01 correlation functionals. 48 The DFT/B3PW91-calculated energy profile agreed well with the CCSD(T)calculated one, as will be discussed in Section 4.1. The (311111/22111/411) basis set was used for Ru with the Stuttgart effective core potentials (ECPs) for replacing its 28 core electrons. 49 Usual 6-311++G(2df,2pd) basis sets were used for O, S, and Cl atoms, 6-311++G(d) basis sets were used for C and N atoms, and 6-31++G(d,p) basis set was used for H atom. In the QM/MM calculation of molecular crystal, we need the LJ parameters. Those of Ru were determined here to be σRu = 2.963 Å and εRu = 0.560 kcal mol−1 by fitting the B3PW91 potential energy surface of ruthenium sulfur dioxide complex with one water (H2 O) molecule. The LJ parameters of H2 O and other atoms were taken from the SPC model 50 and AMBER, 51 respectively. In the switching function of Eq. (5), rinner and router were taken to be 30 and 40 Å, respectively. For making comparison of reaction between solution and crystal phases, we investigated the isomerization of [RuII (NH3 )4 (SO2 )L]n+ (n = 1 for L = Cl− ; n = 2 for L = H2 O (W) and pyridine (Py)) in aqueous phase, employing the 3D-RISM-SCF method with the Kovalenko-Hirata (KH)

10


Page 11 of 55 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


closure. 52–55 As the initial step of the 3D-RISM-SCF-DFT calculation, the 1D-RISM calculation was carried out to set up the solvent-solvent total correlation functions with the density ρW = 1.00 g cm−3 and temperature T = 298.15 K. The point charges and the LJ parameters of H2 O solvent were taken from the SPC model, 50 where the LJ parameters of the water H atom is modified to (σH , εH ) = (1.0 Å, 0.056 kcal mol−1 ). 56–58 The 1D-RISM calculation was performed on a grid of 8096 points, corresponding to the 1D space radius of 56.2 Å. In the 3D-RISM calculation, the LJ parameters of solute molecule were taken to be the same as those employed in the QM/MM calculation of molecular crystal. As a box size of solvent, a cubic grid of 128 points/axis was set with spacing of 0.5 Å for geometry optimization. To divide solvent grid space into three regions, the boundary distances of rin and rout were taken to be 15.0 and 18.0 Å, respectively; see Ref. [32] for the definition. In the RISM-SCF calculation, analytical treatment of the electrostatic asymptotics of the 3D site correlation functions and electrostatic potentials was performed with the nonperiodic asymptotics under the identically zero concentration of ionic species. On such condition, (as)

the long-range asymptotics of the 3D site total correlation function hs (r), can be determined by (as)

evaluating hs (k = 0) with the polynomial extrapolation on a radial k-grid. 59–61

Results and Discussion Thermal Isomerization in Gas Phase Geometries of Ground State and Meta-Stable Isomers and their Relative Stabilities Prior to the discussion of the isomerization in crystal, we investigated the effects of the isomerization of target monomer [RuII (NH3 )4 (SO2 )L]n+ in gas phase, because the isomerization in gas phase is reference to those in aqueous phase and crystal. For convenience, [RuII (NH3 )4 (SO2 )L]n+ with L = Cl− , H2 O (W), and pyridine (Py) is named 1Cl , 1W , and 1Py , respectively, hereafter. To determine the geometry, we employed R-DFT method with the B3PW91 functional because this method provided the energy profile close to that calculated by the CCSD(T) method, as shown in

11



Table S1 in Electronic Supporting Information (ESI); for instance, the DFT/B3PW91-calculated stabilities of MS2 and MS1 relative to GS agree well with those calculated by the CCSD(T). 62 In MS2, two isomers were optimized, as shown in Figure 1 and Scheme 1; in both, the O1 of SO2 takes the position above the Ru atom, where O1 etc. are defined in Scheme 5. One isomer MS2A has the SO2 moiety staggered to the N-Ru-N plane, but the other MS2B has that eclipsed to the N-Ru-N plane. 63 Corresponding to these isomers, two isomers could be optimized for MS2 and GS and two kinds of geometry change could be found here; one is MS1A → TSMM A → MS2A → TSGM A → GSA (path I) and the other is MS1B → TSMM B → MS2B → TSGM B → GSB (path II). 64 The important geometrical parameters determined by the R-DFT/B3PW91 are shown in Figure 1, and the energy profiles were re-calculated by the CCSD(T) in gas phase, as shown in Figure 2A; see Tables S2 and S3 in ESI for details of geometrical parameters and energy profiles. Several characteristic features shown in Figure 1 are summarized, as follows: (i) In MS1, the O1 atom of SO2 interacts with the Ru in the η 1 manner, in which the Ru-O1 distance (2.11 Å in MS1A (path I) and 2.12 Å in MS1B (path II)) is significantly shorter than the Ru-S distance (3.35 Å in MS1A and 3.38 Å in MS1B) by about 1.2 Å. (ii) In TSMM going from MS1 to MS2, the S of SO2 is approaching the Ru by 0.19 Å in path I and 0.31 Å in path II and simultaneously the O1 atom is leaving the Ru by 0.20 Å in path I and 0.16 Å in path II. Also, the Ru-O1 distance of TSMM (2.31 Å in path I and 2.28 Å in path II) is considerably shorter than the Ru-S1 one (3.16 Å in path I and 3.07 Å in path II) by about 0.8 Å. (iii) In MS2, the S=O1 double bond coordinates with the Ru in the η 2 side-on manner, in which the Ru-O1 distance (2.17 Å in MS2A (path I) and 2.18 Å in MS2B (path II)) is somewhat shorter than the Ru-S distance (2.52 Å in MS2A and 2.47 Å in MS2B) by about 0.3 Å. (iv) In TSGM going from MS2 to GS, the Ru-O1 distance (2.86 Å in path I and 2.85 Å in path II) is much longer than that of MS2 but the Ru-S distance (2.50 Å in paths I and II) is almost the same as that of MS2, indicating that the O1 of SO2 is leaving the Ru (by 0.69 Å in path I and 0.67 Å in path II) but the position of the S hardly changes. (v) In GS, only the S atom of SO2 interacts with the Ru in the η 1 manner, in which the Ru-O1 distance (3.15 Å in paths

12


Page 12 of 55

Page 13 of 55 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


I and II) is significantly longer than the Ru-S one (2.18 Å in GSA (path I) and 2.17 Å in GSB (path II).). (vi) The Ru-O1 and Ru-S distances do not monotonously change when going from MS1 to GS but the difference δ r = r(Ru-S) − r(Ru-O1 ) between them monotonously changes (Figure 1 and Figure S1 in ESI), indicating that the difference δ r can be taken as a reaction coordinate of this thermal isomerization. (vii) The above-mentioned characteristic changes of the Ru-O1 and Ru-S distances are similar among 1Cl , 1W , and 1Py (see Table S2 in ESI). (vi) The energy profile largely differs among 1Cl , 1W , and 1Py despite of the similar geometry changes, as shown in Figure 2A; in both of paths I and II, the relative energy ∆Egas increases in the order GS < MS2 ≪ MS1 in 1Cl but GS < MS1 < MS2 in 1W and 1Py . The results of 1W and 1Py indicate that MS1 can be converted to GS by elevating temperature but MS2 cannot be formed from MS1 as stable species, which is completely different from the experimental result in crystal that MS2 was observed at low temperature (e.g. 100 - 140 K in 1Py ) and it is converted to GS by elevating temperature (e.g. 150 - 200 K in 1Py ). (ix) In 1Cl , MS2 is more stable than MS1, indicating that MS2 is formed as a more stable intermediate than MS1 in the thermal isomerization. (x) On the other hand, the ‡ energy barrier Egas from MS1 to MS2 and that from MS2 to GS do not differ very much among

1Cl , 1W , and 1Py ; the former (11 kcal mol−1 in 1Cl , 10 kcal mol−1 in 1W , and 8 kcal mol−1 in 1Py ) moderately increases and the latter (20 kcal mol−1 in 1Cl and 19 kcal mol−1 in 1W and 1Py ) slightly increases in the order 1Py ≤ 1W < 1Cl . To elucidate the reason why the energy profiles of 1W and 1Py are similar to each other but significantly differ from that of 1Cl , we performed the localized molecular orbital energy decomposition analysis (LMO-EDA); 65 see pages S7-S9 and Tables S4 and S5 in ESI for detail. This EDA shows three important features, as follows: (i) In TSGM , TSMM , and MS1, the electrostatic (ES) term ∆EES is responsible for the difference in relative stability among 1Cl , 1W , and 1Py ; In 1Cl , the ES repulsion between the negatively charged Cl− and the O of SO2 is responsible for destabilization in energy at MS1 and TSMM , whereas such ES repulsion of neutral W and Py with the O of SO2 is very small in 1W and 1Py , respectively. (ii) In MS2, the sum of polarization, charge transfer, and mixing terms ∆EPOL + ∆ECT + ∆EMIX is responsible for the difference in stability

13



among 1Cl , 1W , and 1Py . (iii) In TSGM and MS2, the correlation term ∆ECORR also contributes to the difference among them, where ECORR is the difference in energy between Hartree-Fock and DFT calculations, which was named “DFT dispersion energy" in Ref. [65]. In summary, the geometry of the Ru-SO2 moiety depends little on the ligand (L = Cl− , W, and Py) at the position trans to SO2 in gas phase but the stabilities of MS2 and MS1 relative to GS significantly depend on the ligand; the relative energy ∆Egas increases in the order GS < MS2 ≪ MS1 in 1Cl and GS < MS1 < MS2 in 1W and 1Py . NBO Charge Re-Distribution We wish to focus here on the charge re-distribution of ruthenium sulfur dioxide complex when going from GS to MS1, where the NBO analysis was employed. 66 Charge re-distribution occurs in a similar manner among 1Cl , 1W , and 1Py , as shown in Table 1; see also Table S6 in ESI for the detailed change in paths I and II and Table S7 in ESI for those in aqueous phase and crystal. When going from GS to TSGM , the SO2 charge becomes negative by 0.19 - 0.20 e and the Ru dz2 orbital population decreases, indicating that the σ -donation from the SO2 moiety to the Ru dz2 becomes weaker. This is reasonable because the HOMO (lσ in Scheme 2) of SO2 overlaps with Ru dz2 orbital in GS but this overlap becomes smaller by the Ru-S bond elongation in TSGM . As a result, the Ru positive charge increases. 67 However, the charges of four NH3 ligands change little, suggesting that the cis-positioned ligands are not influenced very much by the coordination change of SO2 in comparison to the trans-positioned ligand L. When going from TSGM to MS2, the sum of dxz and dyz orbital populations decreases the most at MS2 (Table S6 in ESI) and the SO2 charge becomes more negative by 0.12 - 0.13 e, because the π -back-donation from the Ru dπ orbital to ∗ becomes stronger in this step; see Scheme 2 for the π ∗ of SO . Consistent with this the SO2 πout 2 out

enhancement of the π -back-donation, the Ru positive charge increases and the σ -donation from NH3 to Ru becomes stronger to increase the positive charges on four NH3 ligands by 0.06 - 0.07 e. When going from MS2 to TSMM , the SO2 charge becomes less negative by 0.10 - 0.13 e, because the π -back-donation from the Ru dπ orbital to the SO2 becomes weaker in this step. Consistent

14


Page 14 of 55

Page 15 of 55 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


with this weakening of the π -back-donation, the σ -donation from NH3 to Ru becomes weaker and the positive charge of four NH3 ligands largely decreases by 0.12 - 0.13 e. 68 When going from TSMM to MS1, the SO2 charge becomes moderately more negative or less positive by 0.02 - 0.05 e unexpectedly despite that the π -back donation in the η 2 side-on coordination is completely broken when going from TSMM to MS1. 69 The CT from the SO2 to the Ru is weaker in MS1 than in TSMM , because the nout is at lower energy than the nin (Scheme 2). Consistent with this weak CT in MS1, the Ru atomic charge becomes moderately more positive in MS1 than in TSMM . In MS1, the Ru positive charge is considerably larger, the SO2 positive charge is smaller, the NH3 positive charge is moderately smaller than in GS, but the charge of L does not change very much from that of GS. Though the SO2 coordinates with the Ru using the lone pair lσ (HOMO) to induce strong CT from the SO2 to Ru in GS and using the nin in MS1, the CT from nin orbital of the SO2 to Ru is much weaker in MS1 than in GS to increase the Ru positive charge in MS1, as described above. Nevertheless the NH3 positive charge decreases in MS1, indicating that the CT from the NH3 to the Ru is weaker in MS1 than in GS. The unexpected weakening of this CT in MS1 can be explained in terms of the position of the O1 atom in MS1. As shown in Figure 1, the negatively charged O1 atom of the SO2 approaches the Ru dx2 −y2 orbital to raise the dx2 −y2 orbital energy and suppress the CT from the NH3 to the Ru dx2 −y2 . In summary, the changes in charge distribution are reasonably understood in terms of the change in the Ru-SO2 coordinate bond.

Thermal Isomerization in Aqueous Phase Geometry and Free Energy Changes In aqueous phase, the important geometrical parameters determined by the 3D-RISM-SCF-R-DFT method 34,35 are shown in Figure S1 and Table S8 in ESI. 64 The free energy profiles were recalculated by the 3D-RISM-SCF-CCSD(T) in aqueous phase, as shown in Figure 2B and Table S9 in ESI. The geometry changes of target molecule in aqueous solution (Table S8 in ESI) occur in a 15



similar manner to those (Table S2 in ESI) in gas phase, as shown by the similar change in reaction coordinate δ r (Figure S1 in ESI), except for several minor differences; the distances between the Ru and neutral ligands (NH3 , SO2 , W, and Py) become moderately shorter in aqueous phase than in gas phase by about 0.03 - 0.05 Å because water solvent enhances polarization of NH3 , SO2 , W, and Py ligands, whereas only the Ru-Cl distance of 1Cl is longer in aqueous phase than in gas phase because the Ru-Cl electrostatic attraction decreases in aqueous phase by screening effect of water solvent. Also, the torsion angle τ (S-O1 -Ru-N) of TSMM A in 1Py is somewhat different between the reactions in gas and aqueous phases, depending on the hydrogen bond formation; the O2 atom of SO2 in TSMM A forms intra-molecular hydrogen bond with one NH3 ligand on the cis-position in gas phase but intermolecular hydrogen bond with solvent molecules in aqueous phase. Despite of the similar geometrical changes by the isomerization, stabilities of MS2 and MS1 relative to GS significantly change in aqueous phase, as shown in Figure 2B and Table 2; in both of paths I and II, the relative free energy ∆A increases in the order GS ≪ MS2 ≪ MS1 in 1Cl , 1W , and 1Py unlike ∆Egas in gas phase (GS < MS2 ≪ MS1 in 1Cl and GS < MS1 < MS2 in 1W and 1Py ). These results mean that, in aqueous phase, MS1 is converted to MS2 and then MS2 is converted to GS in a stepwise manner by increasing temperature. This is different in the gas phase reaction of 1W and 1Py ; remember that MS2 cannot be found as stable species in the isomerization of 1W and 1Py in gas phase because MS2 is less stable than MS1. The free energy barrier A‡ for ‡ TSMM in solution differs little from the Egas in gas phase but the A‡ for TSGM becomes moderately ‡ larger than the Egas by 1 - 2 kcal mol−1 , as shown in Table 2. Because the solute energy Esolute

increases in the order MS1 ≤ GS < MS2 in 1W and 1Py , as shown in Table 3 and Figure S2 in ESI, it is likely concluded that the solvation free energy ∆µ is responsible for increasing order of free energy GS < MS2 < MS1, which is different from the order in gas phase (see also Table S9 in ESI); note that the free energy A is the sum of the solute energy Esolute and the solvation free energy ∆µ . Actually, the ∆µ is about 2 times larger in MS2 and MS1 of 1Py and 1W than that of 1Cl . The same trends of solvation free energy were provided by the polarizable continuum model (PCM) and solvation model density (SMD) model, 71,72 as shown in Table S10 in ESI. The reasons

16


Page 16 of 55

Page 17 of 55 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


why the solvation free energy ∆µ largely differs among GS, MS2, and MS1 and largely depends on the trans-ligand L will be discussed below. In summary, the geometry changes in aqueous phase differ little from those in gas phase but the stabilities of MS2 and MS1 relative to GS are quite different from those in gas phase. The solvation free energy plays crucial role for the relative stabilities of these complexes in aqueous phase.

Analysis of Solvation Effects on Energy Change In aqueous phase, the free energy profiles of 1W and 1Py are similar to that of 1Cl unlike in gas phase, as shown in Figure 2B, though the solute energy profiles of 1W and 1Py (see Figure S2 in ESI) largely differ from that of 1Cl like in gas phase (see Figure 2A). Because the LMO-EDA indicates that the change in the electronic structure differs little between gas and aqueous phases, as shown in pages S7-S9 in ESI, the difference in energy profile between gas and aqueous phases (Table 2) arises from solvation free energy, as follows: MS1 of 1W and 1Py is moderately higher in energy in gas phase than GS by about 4 and 1 kcal mol−1 , respectively, but these differences significantly increase to 20 and 16 kcal mol−1 by solvation effect in aqueous phase, as shown in Figure 2 and Table 2. In 1Cl , on the other hand, the difference in energy between MS1 and GS is calculated to be 15 kcal mol−1 in gas phase and 25 kcal mol−1 in aqueous phase (Figure 2 and Table 2). This means that the solvation free energy ∆µ much more destabilizes MS1 in 1W and 1Py than in 1Cl , as shown in Table 3 and Table S9 in ESI. To correctly understand this reaction in aqueous phase, therefore it is of considerable importance to make analysis of solvation free energy. For the analysis, the linear response solvation free energy Eint /2 of MS1 and its difference (∆Eint /2) from that of GS are evaluated, as shown in Table 2, where Eint /2 is a half of averaged solute-solvent interaction energy. 73 The Eint /2 value of MS2 is not different very much among 1Cl , 1W , and 1Py , as shown in Table 2, but the Eint /2 value of MS1 is much larger in 1W and 1Py than in 1Cl , which is main reason for destabilizing MS1 of 1W and 1Py relative to GS and MS2 in aqueous phase. The Eint is divided into two components, Eint (SO2 ) and Eint (Ru(NH3 )4 L), which are the

17



interactions between SO2 and water solvent and between [Ru(NH3 )4 L]n+ and water solvent (Table 3), respectively. In MS2, the Eint (SO2 )/2 is much more positive in 1W and 1Py than in 1Cl but the Eint (Ru(NH3 )4 L)/2 is significantly more negative in 1W and 1Py than in 1Cl . As a balance of these two terms, the ∆Eint /2 value for MS2 is not very different among 1Cl , 1W , and 1Py . In MS1, the Eint (SO2 )/2 value is similar among these complexes but the Eint (Ru(NH3 )4 L)/2 value is more positive in 1W and 1Py than in 1Cl . These results indicate that the larger destabilization of MS1 in 1W and 1Py than in 1Cl mainly arises from the interaction between the [RuII (NH3 )4 L]n+ moiety and solvent (water) molecules rather than that between the SO2 moiety and solvent molecules. The reasons why the solvation free energy ∆µ is considerably larger (more positive) in MS1 than in GS and MS2 and why the destabilization of MS1 by ∆µ is larger in 1W and 1Py than in 1Cl can be explained by the changes in charge re-distribution of [Ru(NH3 )4 (SO2 )L]n+ and solvation structure when going from GS to MS1; see Figure S9 for the charge re-distribution in aqueous phase, which we will not discuss here because it is similar to that in gas phase (Section 4.1.2). Four important differences in solvation structure were found between GS and MS1, as follows: (i) The distribution of the O site, gO (r), of water solvent is small but the distribution of the H site, gH (r), is large around the O1 atom (blue region a; Scheme 6A left-figure) of SO2 in GS because the O1 atom is negatively charged. In MS1, the gH (r) value largely decreases around the O1 atom of SO2 because the O1 atom coordinates with Ru and the solvation space around the O1 atom becomes very small. (ii) The gO (r) value largely increases around the S atom (red region b) of SO2 in MS1 because the S atom is positively charged and the solvation space around the S atom is small in GS but enough in MS1. (iii) The gH (r) value is moderately smaller around the O2 atom of SO2 (blue region c) and the gO (r) is moderately smaller around two NH3 ligands in MS1 than in GS, because the O2 atom approaches the positively charged H atom of NH3 ligand in MS1 to form the intra-molecular hydrogen bond and hence the solvation space decreases there. (iv) Because the negatively charged O1 atom of SO2 becomes closer to the Ru atom and the positive charge of the NH3 ligand becomes smaller in MS1 than in GS (Table S7 in ESI), 76 the gO (r) value decreases in the red region d (Scheme 6A) around four NH3 ligands. These solvation changes (i),

18


Page 18 of 55

Page 19 of 55 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


(iii), and (iv) lead to weakening of solvation effects to destabilize MS1 relative to GS in aqueous solution, because these solvation changes of (iii) and (iv) decrease the stabilization by the solutesolvent interaction between the [RuII (NH3 )4 ] moiety and water solvent and that of (i) decrease the stabilization between the SO2 moiety and water solvent. Though the solvation change of (ii) leads to enhancement of solvation effects to stabilize MS1 relative to GS, its contribution is not very large. As a result, the stabilization by solvation is smaller in MS1 than in GS in aqueous phase. It is noted here that the difference in ∆µ between MS1 and GS is, however, intrinsically small in 1Cl , because the gH (r) value increases in the region d to destabilize the positively charged [RuII (NH3 )4 ]n+ moiety but to stabilize the negatively charged Cl− . In 1W and 1Py , on the other hand, such cancellation does not occur because neutral W and Py ligands exist instead of negatively charge Cl− . Therefore, the solvation free energy becomes much more positive (destabilizing) in MS1 of 1W and 1Py than in that of 1Cl . It is likely that the above-mentioned difference in solvation free energy between GS and MS1 is general phenomenon in transition metal complex. For instance, the similar solvation is found in the interaction between water molecule (W) and square planar Pt complex [PtII (NH3 )4 ]2+ , which was recently investigated with 3D-RISM-SCF-MP2 method. 77 In H2 O interaction with metal center, there are two possible structures; in one, the O atom of water approaches the metal as expectedly, and in the other, the H atom of water approaches the metal unexpectedly. The former is called O-interaction form and the later is H-interaction form. From the viewpoint of the charge redistribution, the H-interaction and O-interaction forms correspond to GS and MS1, respectively, as shown in Scheme 6B, because in the former the positively charged H approaches Pt and in the latter the negatively charged O approaches Pt. In the case of Pt complex, the O-interaction form of W is less stabilized by water solvent than the H-interaction form like MS1, where W exists on the z-axis. 77 The difference in solvation between the H- and O-interaction forms can be explained in the same way as here. 78 The other factor having influence on ∆µ is the CT from NH3 ligand to Ru, which is weaker in MS1 than in GS, as was discussed above (Table 1). The S atom of SO2 in [RuII (NH3 )4 (SO2 )L]n+

19



is positively charged and the O atom is negatively charged. In GS, the S atom approaches the Ru to lower the Ru dx2 −y2 orbital energy by electrostatic interaction, which enhances CT from NH3 ligand to Ru, as schematically shown in Scheme 6B. In MS1, on the other hand, the O atom approaches the Ru to raise the Ru dx2 −y2 orbital energy, which suppresses the CT from NH3 ligand to Ru (Scheme 6B). As a result, positive charge of NH3 decreases, which weakens the solutesolvent ES interaction between water solvent and four NH3 ligands in MS1 compared to that in GS. Because of these differences in solvation structure between MS1 and GS, MS1 becomes less stable than GS in aqueous phase. The destabilization of MS1 by solvation effect is not very large for anionic L (Cl− ) but large for neutral W and Py, because [RuII (NH3 )4 (SO2 )Cl]+ has one plus charge but [RuII (NH3 )4 (SO2 )Py]2+ and [RuII (NH3 )4 (SO2 )W]2+ have two plus charges; remember that the solvation effect on species with one positive charge is smaller than that on species with two positive charges. In aqueous phase, therefore, the free energy profile of the isomerization reaction of 1Cl changes less by solvation than those of 1W and 1Py .

Thermal Isomerization in Crystal Geometry changes and Energy Profiles In the theoretical study of the thermal isomerization in crystal, we employed [Ru(NH3 )4 (SO2 )Py]2+ (2Cl− ·H2 O) because the large quantum yield of photo-reaction was reported and meta-stable states were experimentally observed well at higher temperature (around 100 - 140 K) in this complex than in other ruthenium sulfur dioxide complexes; 29 remember that the use of temperature close to cryogenic temperature is not good for application. The potential energy profile of [Ru(NH3 )4 (SO2 )Py]2+ (2Cl− ·H2 O), where 2Cl− are counterion and H2 O is crystal water, was calculated by the QM/MM method, as shown in Figure 3; see brackets in Figure S1 in ESI for the geometry change in the MS1-crystal at 140 K and Table S11 in ESI for details of geometry in other crystal models. 64 As shown by the similar values of reaction coordinate δ r among gas, aqueous, and crystal 20


Page 20 of 55

Page 21 of 55 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


phases (Figure S1 in ESI), the geometry changes (Table S11 in ESI) of the target molecule 1Py in crystal occur in a similar manner to those along path I in gas and solution phases (Tables S2 and S4 in ESI), except for three characteristic differences; (i) geometry changes corresponding to path II disappears because of large steric repulsion between SO2 moiety and counter Cl− anion on the x-axis. (ii) The Ru-O1 distance of MS1 is much shorter in crystal phase than in aqueous phase by about 0.05 - 0.08 Å, which is induced by the ES and steric repulsions between the SO2 of the target QM molecule and the neighboring Ru complex on the +z-axis (Scheme 7C), as will be discussed in Section 4.3.2. (iii) The torsion angle τ (S-O1 -Ru-N) of TSMM in the MS2-crystal is somewhat different from that in the MS1- and GS-crystals; the τ (S-O1 -Ru-N) of TSMM in the MS1-crystal is almost the same as that of TSMM A in aqueous phase in which the O2 atom of SO2 forms the intermolecular hydrogen bond with surrounding molecules. On the other hand, the τ (S-O1 -RuN) of TSMM in the MS2-crystal is closer to that of TSMM A in gas phase than in aqueous phase, indicating that TSMM in the MS2-crystal has the tendency to form the intra-molecular hydrogen bond between the SO2 and one NH3 ligand on the cis-position. Despite of the similar geometrical changes by the isomerization, stabilities of MS2 and MS1 relative to GS significantly change in crystal phase compared to those in gas phase; in crystal, MS2 exists at much higher energy than GS and MS1 exists at much higher energy than MS2 unlike in gas phase, independently of the lattice vector and the kinds of crystal model, as shown in Figure 3 and Table S12 in ESI. This trend is similar to the free energy profile in aqueous phase (Figure 2B). However, the energy barrier of the MS2→GS isomerization and that of the MS1→MS2 isomerization depend on both of the lattice vector and the kinds of crystal model, as shown by E ‡ values in three kinds of crystal in Figure 3. MS2 is the most stable relative to GS in the MS2and MS2′ -crystals and the least stable in the GS-crystal. 79 MS1 is the most stable in the MS1crystal but less stable in the GS- and MS2-crystals. These results about the stability dependence on crystal structure is consistent with our intuitive understanding, to which theoretical support has not been reported. We wish to emphasize that the above-mentioned computational result is the first theoretical evidence showing that molecule is the most stable in the crystal consisting of the same

21



Page 22 of 55

molecules. To save space, we omit here the discussion about the thermal isomerization in the MS2′ -crystal, because the potential energy profile differs little between the MS2- and MS2′ -crystals (Figure 3). Hereafter, we will discuss the MS1→MS2 isomerization in the MS1-crystal and the MS2→GS isomerization in the MS2-crystal, because it is likely that most of MS1 molecules are surrounded by MS1-crystal and most of MS2 molecules are surrounded by MS2-crystal when the isomerization starts to occur. 29 Figure 3b indicates that the energy barrier for the MS2→GS isomerization depends little on the lattice vector but that for the MS1→MS2 isomerization somewhat decreases in the lattice vector at higher temperature (140 K). 46 Judging from the energy difference between MS1 and TSMM , the MS1→MS2 isomerization occurs with activation barrier E ‡ of 10.9 kcal mol −1

in the lattice vector at 100 K and 8.6 kcal mol

−1

in the lattice vector at 140 K (Figure 3b).

The MS2→GS isomerization needs the barrier of 17.1 mol

−1

in the lattice vector at 100 K and

140 K (Figure 3c). It is noted that the E ‡ of the MS1→MS2 isomerization is considerably smaller than that of the MS2→GS isomerization. These results lead to the conclusion that the MS1→MS2 isomerization occurs at lower temperature than the MS2→GS isomerization, which is consistent with the experimental observation that the MS1→MS2 isomerization occurs at low temperature 100 - 140 K but the MS2→GS isomerization does not occur at the same temperature but needs higher temperature 150 - 200 K. 29

Analysis of Crystal Effects on Energy Profiles and Characteristic Features of Reaction in Crystal The relative stabilities of GS, MS2, and MS1 of 1Py , and the activation energy of its isomerization are compared among gas, aqueous, and crystal phases in Table 4; see also Table S13 in SI for the rate constant. The order of rate constant for the isomerization in crystal agrees with the experimental observation that the MS1→MS2 and MS2→GS isomerizations occur completely at about 140 and 200 K. respectively. 29 Interestingly, the activation barriers in crystal are similar to those in gas phase and aqueous phase, but the relative stabilities of meta-stable isomers (MS1 and

22


Page 23 of 55 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


MS2) considerably differ among gas, aqueous, and crystal phases and therefore the energetics for isomerization considerably differs among those there phases. We investigated the effects of the crystal by separating the surrounding crystalline molecules into the nearest neighboring QM1 (two neighboring Ru complexes 1Py ) and QM2 (five Cl− and one H2 O), the next neighboring eight ruthenium complexes (MM1 to MM8), and the long-range ES interaction with remaining MM crystal, as shown in Scheme 4. The influences of QM1 and QM2 on the target monomer QM0 are shown in Table 5 and Figure S3 in ESI; In this Table, 0+1, 0+2, and 0+1+2 represent QM systems consisting of QM0 and QM1, QM0 and QM2, and QM0, QM1, and QM2, respectively, where the geometries of QM0, QM1, and QM2 are taken to be the same as those optimized by the QM/MM method in Figure 3. In 0+1+2, the QM calculation was performed in the presence of all MM molecules in MS1-crystal. 80 As shown in Table S14 in ESI, the potential energy E(QM0) changes in almost the same manner among GS-, MS1-, and MS2-crystals when going from MS1 to GS, which is similar to that in gas phase (Figure S1 in ESI); for instance, E(QM0) of MS1 is somewhat below that of MS2 by about 3 kcal mol−1 in all kinds of crystal model and the differences in E(QM0) between MS1 and TSMM and between MS2 and TSGM depend little on the lattice vectors at 100, 120, and 140 K. These results lead to the conclusion that crystal effects from QM1, QM2, and MM crystal are crucially important to understand correctly reaction in crystal. The effects of QM1, QM2, MM1-MM8, and remaining part of MM crystal are compared in MS1→MS2 isomerization in the MS1-crystal and the MS2→GS isomerization in the MS2crystal, as shown in Table 5. Those results are summarized, as follows: (i) The sum of these crystalline effects is much larger (about 2 times larger) than the solvation effects. (ii) The effects of QM1 and QM2 are not very large for GS and TSGM . For instance, they destabilize TSGM relative to GS by 2 - 6 kcal mol−1 , which is smaller than that by the solvation effect (5 - 6 kcal mol−1 ), as shown by the values of 0+1 and 0+2 in Table 5; see also Figures 2 and S3 in ESI. (iii) QM2 largely destabilizes MS2 by 9 - 12 kcal mol−1 relative to GS, whereas QM1 moderately destabilizes MS2 by 2 - 3 kcal mol−1 relative to GS, which is similar to that (3 - 4 kcal mol−1 )

23



by the solvation effect. This large effect of QM2 is reasonable because the O2 atom of the SO2 moiety is closer to one of Cl− counterions in MS2 than in GS and MS1 (remember the η 2 side-on coordinated structure of SO2 in MS2) and therefore the ES repulsion between the O2 atom of SO2 and Cl− counterion is the largest in MS2, as shown in Schemes 7A and 7B. These are also the reasons why MS2 exists at higher energy in crystal than in gas and aqueous phases (Table 4). (iv) The effects of QM1 and QM2 are very large for TSMM and MS1; for instance, QM1 significantly destabilizes TSMM and MS1 by 9 - 10 kcal mol−1 and 13 - 16 kcal mol−1 , respectively, and QM2 significantly destabilizes them by 13 - 15 kcal mol−1 and 14 - 15 kcal mol−1 , respectively. It is noted that the sum of these effects is almost 2 times larger than those by the solvation effect (about 15 kcal mol−1 in TSMM and 14 kcal mol−1 in MS1). Such considerably large energy destabilizations of TSMM and MS1 by QM1 is reasonable because the neighboring Ru complex 1Py on the +z-axis (QM1) has the positive charge and the ES attraction between this positively charged 1Py (QM1) and the O atom of SO2 (QM0) becomes weaker in MS1 and TSMM relative to GS; remember that both of O1 and O2 atoms of the SO2 (QM0) are close to the neighboring 1Py on the +z-axis (QM1) in GS but distant in MS2 and MS1 (Scheme 7). In addition to the ES effect of QM1, the distance between the S of SO2 in QM0 and HPy of QM1 is short to induce the steric repulsion and destabilize TSMM and MS1. The S· · · HPy distance (2.17 Å) in the MS1-crystal is considerably shorter for the lattice vector at 140 K than that (2.33 Å) for the lattice vector at 100 K (see Scheme 7 and Table S11 in ESI), indicating that the steric repulsion between QM1 and the SO2 of QM0 is responsible for the smaller energy barrier for TSMM in the lattice vector at 140 K that that at 100 K. QM2 (five Cl− and one H2 O) also contributes to the crystalline effect on TSMM and MS1 because the positively charged S atom of SO2 (QM0) is close to two Cl− (QM2) in GS but the negatively charged O1 atom of the SO2 (QM0) is close to them in MS1 (Schemes 7A and 7C). (v) MM1 to MM8 and the remaining part of the crystal do not influence very much the energetics compared to QM1 and QM2 (Table 5), because most of them are canceled by each other, as shown in Figure S13 in ESI. This means that the effects of QM1 and QM2 are mainly responsible for the relative stabilities of GS, MS2, and MS1 in the crystal. However, the important

24


Page 24 of 55

Page 25 of 55 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


roles of periodically infinite MM crystal is clearly shown by model calculations with the absence of the infinite MM crystal. 80 Also, the presence of MM1 to MM8 in addition to QM1 and QM2 is not enough to reproduce the energy barriers for TSGM and TSMM , especially, in the lattice vector at 140 K, as shown in Table 5; see Figure S4 in ESI for more details. On the basis of these results, it is concluded that the incorporation of the short-range and long-range ES interactions by charged species is crucially important to describe correctly characteristic feature of molecule in crystal. In summary, the nearest neighboring molecules and counter anions play crucially important roles in crystalline effects and the periodically infinite MM crystal is also important to understand correctly the reaction in crystal. In the solution reaction, solvent position and orientation can change easily, according to the solute geometry changes. Though the solvation effects are enough strong to change the stability order of GS, MS1, and MS2, the solvation effects are smaller than the crystalline effects in which crystal geometry is not flexible expectedly.

Conclusions Though DFT calculations with plane wave basis set have been often employed under periodic boundary conditions, the use of hybrid DFT functional and post Hartree-Fock method is still difficult for theoretical study of crystal even nowadays. In particular, it is very difficult to perform QM study of chemical reaction in solid crystal. We investigated the surrounding molecular effects on the thermal isomerization of ruthenium sulfur dioxide complex [RuII (NH3 )4 (SO2 )(pyridine)]2+ 1Py in crystal, using our QM/MM method with DFT and CCSD(T) methods. To clearly show characteristic features of reaction in crystal, we compared geometry and energy changes by the isomerization reaction of this ruthenium complex in crystal with those in gas and solution phases. The isomerization thermally occurs from the η 1 -O-bound SO2 complex MS1 to the η 1 -S-bound SO2 complex GS through η 2 -SO-bound SO2 complex MS2. We optimized transition state TSMM connecting MS1 and MS2 and TSGM connecting MS2 and GS in gas, water, and crystal phases. In gas phase, [Ru(NH3 )4 (SO2 )Cl]+ 1Cl becomes lower in energy following the order MS1 >

25



MS2 > GS. However, [Ru(NH3 )4 (SO2 )(H2 O)]2+ 1W and 1Py become lower in energy following the order MS2 > MS1 > GS, indicating that MS2 cannot be observed because the MS2→MS1 isomerization occurs more easily than the MS1→MS2 isomerization in 1W and 1Py . In aqueous solution, all these complexes become lower in free energy A following the order MS1 > MS2 > GS and the MS1→MS2 isomerization occurs with smaller activation barrier than the MS2→GS one. It is noted that MS1 of 1W and 1Py becomes more unstable than MS2 in aqueous solution. The relative energy of the solute differs little between gas and aqueous phases, indicating that the solvation free energy induces the large energy destabilization of MS1 relative to MS2 and GS. The significantly large destabilization of solvation free energy at MS1 in 1W and 1Py arises from the decrease in the solute-solvent ES attraction between the [Ru(NH3 )4 L]2+ moiety and solvent (water) molecules in MS1. The decrease in ES attraction arises from the approach of negatively charged O1 atom of SO2 to the Ru, which induces the weakening of CT from four NH3 to the Ru at MS1, the decrease in the gO (r) value around the [Ru(NH3 )4 L]2+ moiety, and the large destabilization of solvation free energy in MS1. We focused on the thermal isomerization of 1Py in molecular crystal because this system exhibits the large quantum yield for photo-isomerization and its thermal isomerization has been experimentally investigated well. We successfully optimized the important five geometries, GS, TSGM , MS2, TSMM , and MS1, incorporating the crystalline effects of GS-, MS1-, MS2-, and MS2′ -crystals using our periodic QM/MM method. In all crystal models, 1Py becomes lower in energy following the order MS1 > MS2 > GS unlike in gas phase. Both stabilities of MS1 and MS2 relative to GS moderately depend on the kinds of the crystal structure and lattice vector. The MS1→MS2 isomerization occurs with smaller activation barrier than the MS2→MS1 one. These results agree with the experimental observation that the decay of MS1 into MS2 occurs at low temperature but the decay of MS2 into GS occur at elevated temperature. To investigate the reason why the energy profile largely differs between gas phase and molecular crystals, we analyzed the influences of surrounding molecules on the target QM molecule. The geometry changes of the target molecule in crystal occur in a similar manner to those in gas and

26


Page 26 of 55

Page 27 of 55 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


solution phases, indicating that insight only into the target molecule is not meaningful at all and the crystalline effects are crucially important for understanding reaction in crystal. Two neighboring Ru complexes on the z-axis largely destabilize MS1 relative to GS and five neighboring Cl− counterions largely destabilize MS2 and MS1 relative to GS. These results are reasonable because the nearest neighboring Ru complex (QM1) on the +z-axis induces the ES repulsion between the positive net charge of QM1 and the positively charged S atom of the SO2 in QM0 and also five neighboring Cl− counterions (QM2) induce the ES repulsion with the negatively charged O atom of the SO2 in QM0. In addition, the steric repulsion between the SO2 moiety of QM0 and the pyridine moiety of the nearest neighboring Ru complex is responsible for the large destabilization of MS1 and small energy barrier of the MS1→MS2 isomerization in the MS1-crystal. However, the effects by these two neighboring Ru complexes on the z-axis and five neighboring Cl− counterions are not enough to explain quantitatively the energy change by the isomerization. The long-range ES interaction between the target QM0 and periodically infinite MM crystal must be taken into consideration to correctly evaluate the energy profile in crystal. Therefore, the present QM/MM method based on the periodic boundary condition is useful to investigate chemical reaction in molecular crystal. In the present study, our QM/MM method was successfully applied to the thermal isomerization in crystal, as discussed above. Our QM/MM method can be applied to unimolecular photoisomerization in crystal using CASSCF/CASPT2, EOM-CC, and TD-DFT methods for QM calculation.

Acknowledgement This work was financially supported by Grants-in-Aid for Young Scientists (No. 17K14437) and for Scientific Research (No. 18H04513, and No. 15H03770). This paper is based on some results obtained from a project commissioned by the New Energy and Industrial Technology Development Organization (NEDO) supported by the Ministry of Economy, Trade and Industry, Japan. We are grateful to IMS computer center and ACCMS, Kyoto University for the supply of CPU time. 27



Supporting Information Available Details of RISM-SCF method, complete reference of [42], energy decomposition analysis, dependence of energy profile on computational methods, important geometrical parameters, change in NBO charge, illustration of MM1 to MM8, and influences of surrounding QM1, QM2, and MM1 to MM8 on target Ru complex QM0. This material is available free of charge via the Internet at http://pubs.acs.org/.

28


Page 28 of 55

Page 29 of 55 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


References (1) Coppens, P.; Novozhilova, I.; Kovalevsky, A. Photoinduced linkage isomers of transitionmetal nitrosyl compounds and related complexes. Chem. Rev. 2002, 102, 861-884. (2) Schaniel, D.; Imlau, M.; Weisemoeller, T.; Woike, T.; Krämer, K. W.; Güdel, H. U. Photoinduced nitrosyl linkage isomers uncover a variety of unconventional photorefractive media. Adv. Mater. 2007, 19, 723-726. (3) Bitterwolf, T. E. Photochemical nitrosyl linkage isomerism/metastable states. Coord. Chem. Rev., 2006, 250, 1196-1207. (4) Coppens, P.; Fomitchev, D. V.; Carducci, M. D.; Culp, K. Crystallography of molecular excited states. Transition-metal nitrosyl complexes and the study of transient species. J. Chem. Soc. Dalton Trans. 1998, 6, 865-872. (5) Imlau, M.; Haussühl, S.; Woike, T.; Schieder, R.; Angelov, V.; Rupp, R. A.; Schwarz, K. Holographic recording by excitation of metastable electronic states in Na2 [Fe(CN)5 (NO)]· 2H2 O: A new photorefractive effect. Appl. Phys. B 1999, 68, 877-885. (6) Cole, J. M. Photocrystallography. Acta Crystallogr. A 2008, 64, 259-271. (7) Kovalevsky, A. Y.; King, G.; Bagley, K. A.; Coppens, P. Photoinduced oxygen transfer and double-linkage isomerism in a cis-(NO)(NO2 ) transition-metal complex by photocrystallography, FT-IR spectroscopy and DFT calculations. Chem. Eur. J. 2005, 11, 7254-7264. (8) Schaniel, D.; Cormary, B.; Malfant, I.; Valade, L.; Woike, T.; Delley, B.; Krämer, K. W.; Güdel, H. U. Photogeneration of two metastable NO linkage isomers with high populations of up to 76% in trans-[RuCl(py)4 (NO)][PF6 ]2 · 12 H2 O. Phys. Chem. Chem. Phys. 2007, 9, 3717-3724. (9) Schaniel, D.; Woike, T.; Delley, B.; Biner, D.; Krämer, K. W.; Güdel, H. U. Generation of one

29



light-induced metastable nitrosyl linkage isomer in [Pt(NH3 )4 Cl(NO)]Cl2 in the red spectral range. Phys. Chem. Chem. Phys. 2007, 9, 5149-5157. (10) Schaniel, D.; Woike, T.; Kushch, L.; Yagubskii, E. Photoinduced nitrosyl linkage isomers in complexes based on the photochromic cation [RuNO(NH3 )5 ]3+ with the paramagnetic anion [Cr(CN)6 ]3− and the diamagnetic anions [Co(CN)6 ]3− and [ZrF6 ]2− . Chem. Phys. 2007, 340, 211-216. (11) Cormary, B.; Malfant, I.; Cointe, M. B.; Toupet, L.; Delley, B.; Schaniel, D.; Mockus, N.; Woike, T.; Fejfarová, K.; Petricek, V.; Dusek, M. [Ru(py)4 Cl(NO)](PF6 )2 · 0.5 H2 O: A model system for structural determination and ab initio calculations of photo-induced linkage NO isomers. Acta Crystallogr. B 2009, 65, 612-623. (12) Schaniel, D.; Mockus, N.; Woike, T.; Klein, A.; Sheptyakov, D.; Todorova, T.; Delley, B. Reversible photoswitching between nitrito-N and nitrito-O isomers in trans-[Ru(py)4 (NO2 )2 ]. Phys. Chem. Chem. Phys. 2010, 12, 6171-6178. (13) Goulkov, M.; Schaniel, D.; Woike, T. Pulse recording of thermal and linkage isomer gratings in nitrosyl compounds. J. Opt. Soc. Am. B 2010, 27, 927-932. (14) Warren, M. R.; Brayshaw, S. K.; Johnson, A. L.; Schiffers, S.; Raithby, P. R.; Easun, T. L.; George, M. W.; Warren, J. E.; Teat, S. J. Reversible 100% linkage isomerization in a single-crystal to single-crystal transformation: Photocrystallographic identification of the metastable [Ni(dppe)(η 1 -ONO)Cl] isomer. Angew. Chem. Int. Edn 2009, 48, 5711-5714. (15) Hatcher, L. E.; Warren, M. R.; Allan, D. R.; Brayshaw, S. K.; Johnson, A. L.; Fuertes, S.; Schiffers, S.; Setevenson, A. J.; Teat, S. J.; Woodall, C. H.; Raithby, P. R. Metastable linkage isomerism in [Ni(Et4 dien)(NO2 )2 ]: A combined thermal and photocrystallographic structural investigation of a nitro/nitrito interconversion. Angew. Chem. Int. Edn 2011, 50, 8371-8374. (16) Brayshaw, S. K.; Easun, T. L.; George, M. W.; Griffin, A. M. E.; Johnson, A. L.; Raithby, P. R.; Savarese, T. L.; Schiffers, S.; Warren, J. E.; Warren, M. R.; Teat, S. J. Photocrystallo30


Page 30 of 55

Page 31 of 55 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


graphic identification of metastable nitrito linkage isomers in a series of nickel(II) complexes. Dalton Trans. 2012, 41, 90-97. (17) Fomitchev, D. V.; Bagley, K. A.; Coppens, P. The first crystallographic evidence for side-on coordination of N2 to a single metal center in a photoinduced metastable state. J. Am. Chem. Soc. 2000, 122, 532-533. (18) Fomitchev, D. V.; Novozhilova, I.; Coppens, P. Photo-induced linkage isomerism of transition metal nitrosyl and dinitrogen complexes studied by photocrystallographic techniques. Tetrahedron 2000, 56, 6813-6820. (19) Rack, J. J.; Winkler, J. R.; Gray, H. B. Phototriggered Ru(II)-dimethylsulfoxide linkage isomerization in crystals and films. J. Am. Chem. Soc. 2001, 123, 2432-2433. (20) Rachford, A. A.; Petersen, J. L.; Rack, J. J. Designing molecular bistability in ruthenium dimethyl sulfoxide complexes. Inorg. Chem. 2005, 44, 8065-8075. (21) Mockus, N. V.; Marquard, S.; Rack, J. J. Excited state isomerization in a new ruthenium chelating sulfoxide complex. J. Photochem. Photobiol. A 2008, 200, 39-43. (22) Mockus, N. V.; Rabinovich, D.; Petersen, J. L.; Rack, J. J. Femtosecond isomerization in a photochromic molecular switch. Angew. Chem., Int. Ed. 2008, 47, 1458-1461. (23) McClure, B. A.; Abrams, E. R.; Rack, J. J. Excited state distortion in photochromic ruthenium sulfoxide complexes. J. Am. Chem. Soc. 2010, 132, 5428-5436. (24) Kovalevsky, A. Y.; Bagley, K. A.; Coppens, P. The first photocrystallographic evidence for light-induced metastable linkage isomers of ruthenium sulfur dioxide complexes. J. Am. Chem. Soc. 2002, 124, 9241-9248. (25) Kovalevsky, A. Y.; Bagley, K. A.; Cole, J. M.; Coppens, P. Light-induced metastable linkage isomers of ruthenium sulfur dioxide complexes. Inorg. Chem. 2003, 42, 140-147.

31



(26) Bowes K. F.; Cole, J. M.; Husheer, S. L. G.; Raithby, P. R.; Savarese, T. L.; Sparkes, H. A.; Teat, S. J.; Warren, J. E. Photocrystallographic structure determination of a new geometric isomer of [Ru(NH3 )4 (H2 O)(η 1 -OSO)][MeC6 H4 SO3 ]2 . Chem. Commun. 2006, 23, 2448-2450. (27) Phillips, A. E.; Cole, J. M.; d’Almeida, T.; Low, K. S. Effects of the reaction cavity on metastable optical excitation in ruthenium-sulfur dioxide complexes. Phys. Rev. B 2010, 82, 155118. (28) Phillips, A. E.; Cole, J. M.; d’Almeida, T.; Low, K. S. Ru-OSO coordination photogenerated at 100 K in tetraammineaqua(sulfur dioxide)ruthenium(II) (±)-camphorsulfonate. Inorg. Chem. 2012. 51, 1204-1206. (29) Sylvester, S. O.; Cole, J. M.; Waddell, P. G. Photoconversion bonding mechanism in ruthenium sulfur dioxide linkage photoisomers revealed by in situ diffraction. J. Am. Chem. Soc. 2012, 134, 11860-11863. (30) Sylvester, S. O.; Cole, J. M.; Waddelll, P. G.; Nowell, H.; Wilson, C. SO2 phototriggered crystalline nanomechanical transduction of aromatic rotors in tosylates: Rationalization via photocrystallography of [Ru(NH3 )4SO2 X]tosylate2 (X = pyridine, 3-Cl-pyridine, 4-Cl-pyridine). J. Phys. Chem. C 2014, 118, 16003-16010. (31) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227-249. (32) Aono, S.; Sakaki, S. Proposal of new QM/MM approach for geometry optimization of periodic molecular crystal: Self-consistent point charge representation for crystalline effect on target QM molecule. Chem. Phys. Lett. 2012, 544, 77-82. See also pages S2-3 in ESI for the details of RISM-SCF method.

32


Page 32 of 55

Page 33 of 55 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


(33) Yoshida, N.; Hirata, F. A new method to determine electrostatic potential around a macromolecule in solution from molecular wave functions. J. Comput. Chem. 2006, 27, 453-462. (34) Aono, S.; Sakaki, S. Evaluation procedure of electrostatic potential in 3D-RISM-SCF method and its application to hydrolyses of cis- and transplatin complexes. J. Phys. Chem. B 2012, 116, 13045-13062. (35) Aono, S.; Hosoya, T.; Sakaki, S. A 3D-RISM-SCF method with dual solvent boxes for a highly polarized system: application to 1,6-anhydrosugar formation reaction of phenyl α and β -D-glucosides under basic conditions. Phys. Chem. Chem. Phys. 2013, 17, 6368-6381. (36) Iwano, K.; Shimoi, Y. Large electric-potential bias in an EDO-TTF tetramer as a major mechanism of charge ordering observed in its PF6 salt: A density functional theory study. Phys. Rev. B 2008, 77, 075120. (37) Iwano, K.; Shimoi, Y. Revealing the photorelaxation mechanism in a molecular solid using density-functional theory. Phys. Rev. Lett. 2013, 110, 116401. (38) In Ref. [37], geometry optimization with self-consistent point-charges was carried out using standard ONIOM method. In the standard ONIOM method, the calculation at high-level layer (QM) does not explicitly influence that at low-level layer (MM). This condition is not satisfied when the self-consistent point-charges are introduced, because the self-consistent point-charges at the low-level layer (MM) are determined by the calculation at the high-level layer (QM) in an SCF manner. In this regard, the derivatives of self-consistent point-charges with respect to molecular coordinates must be involved in the energy gradient. However, the dependence of self-consistent point-charge on molecular coordinates was not considered in their geometry optimization but it was considered in our method. (39) Remember that Eq. (2) is based on the assumption that all the symmetrically minimal subunits takes the same geometry and electronic structure.

33



(40) In addition to the optimization of transition state, geometries of GS, MS1, and MS2 were refined by the QM/MM calculation under each MM model of GS-, MS1-, and MS2-crystals to compare their relative stabilities at the same computational quality. (41) In our QM/MM calculation of isomerization in crystal, the relaxation of molecules surrounding QM0 would be underestimated in the transition state and product in comparison to the reactant, because the MM part constructed by modeling the crystal consisting of the reactant (Scheme 3a) was used in calculations of transition state and product. Thus, our QM/MM calculation would lead to overestimation of activation energy. However, it is likely that the overestimation is not large because the positions and orientations of the nearest neighboring QM1 and QM2 were optimized in the calculation. (42) Svensson, M.; Humbel, S.; Froese, R. D. J.; Matsubara, T.; Sieber, S.; Morokuma, K. ONIOM: A multilayered integrated MO + MM method for geometry optimizations and single point energy predictions. A test for Diels-Alder reactions and Pt(P(t-Bu)3 )2 + H2 oxidative addition. J. Phys. Chem. 1996, 100, 19357-19363. (43) Dapprich, S.; Komaromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J. A new ONIOM implementation in Gaussian 98. 1. The calculation of energies, gradients and vibrational frequencies and electric Field derivatives. J. Mol. Struct. THEOCHEM, 1999, 462, 1-21. (44) Schmidt, M. W.; Baldridge, K. K.; Boatz J. A. et al. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347-1363. See page S6 in ESI for the complete reference. (45) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran 77, 2nd ed. Cambridge University Press, NewYork, 1992. (46) The neighbor molecules were taken as the QM molecules to improve the description of intermolecular interactions with the target key molecule. In this study, we determined the crystal structure and discussed the energy profile in crystal based on the potential energy without 34


Page 34 of 55

Page 35 of 55 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 entropy contribution. Also, the dependence of reactivity on the lattice vector was taken into consideration as the potential energy changes using the experimentally observed lattice vector at given temperatures 100, 120, and 140 K. (47) Becke, A. D. Density‐functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648-5652. (48) Perdew, J. P.; Wang, W. R. Accurate and simple analytic representation of the electron-gas correlation-energy. Phys. Rev. B 1992, 45, 13244-13249. (49) Andrae, D.; Haeussermann, U.; Dolg, M.; Stoll, H.; Preuss, H. Energy-adjustedab initio pseudopotentials for the second and third row transition elements. Theor. Chem. Acc. 1990, 77, 123-141. (50) Berendsen, H. J. C.; Postma, J. P. M.; van Gunstern, W. F.; Hermans, J. Intermolecular Forces; Pullman, B. Ed.; Reidel, Dordrecht, 1981, pp 331-342. (51) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 1995, 117, 5179-5197. (52) Kovalenko, A. Three-dimensional RISM theory for molecular liquids and solid-liquid interfaces, In Molecular Theory of Solvation; Hirata, F., Ed.; Understanding Chemical Reactivity (series); Mezey, Paul G., Series Ed.; Kluwer Acadamic Publishers: Dordrecht, The Netherlands, 2003; Vol. 24, pp 169-275. (53) Kovalenko, A.; Hirata, F. Three-dimensional density profiles of water in contact with a solute of arbitrary shape: A RISM approach. Chem. Phys. Lett. 1998, 290, 237-244. (54) Kovalenko, A.; Hirata, F. Potentials of mean force of simple ions in ambient aqueous solution.

35



I. Three-dimensional reference interaction site model approach. J. Chem. Phys. 2000, 112, 10391-10402. (55) Kovalenko, A.; Hirata, F. Potentials of mean force of simple ions in ambient aqueous solution. II. Solvation structure from the three-dimensional reference interaction site model approach, and comparison with simulations. J. Chem. Phys. 2000, 112, 10403-10417. (56) Ten-no, S.; Hirata, F.; Kato, S. A hybrid approach for the solvent effect on the electronic structure of a solute based on the RISM and Hartree-Fock equations. Chem. Phys. Lett. 1993, 214, 391-396. (57) Ten-no, S.; Hirata, F.; Kato, S. Reference interaction site model self-consistent field study for solvation effect on carbonyl compounds in aqueous solution. J. Chem. Phys. 1994, 100, 7443-7453. (58) Sato, H.; Hirata, F.; Kato, S. Analytical energy gradient for the reference interaction site model multiconfigurational selfconsistent-field model. J. Chem. Phys. 1996, 105, 1546-1551. (59) Kaminski, J. W.; Gusarov, S.; Wesolowski, T. A.; Kovalenko, A. Modeling solvatochromic shifts using the orbital-free embedding potential at statistically mechanically averaged solvent density. J. Phys. Chem. A 2010, 114, 6082-6096. (60) Gusarov, S.; Pujari, B. S.; Kovalenko, A. Efficient treatment of solvation shells in 3D molecular theory of solvation. J. Comput. Chem. 2012, 33, 1478-1494. (61) Perkyns, J. S.; Lynch, G. C.; Howard, J. J.; Pettitt, B. M. Protein solvation from theory and simulation: Exact treatment of Coulomb interactions in three-dimensional theories. J. Chem. Phys. 2010, 132, 064106. (62) Because the DFT/B3PW91 method reproduced better the CCSD(T)-calculated energetics than the B3PW91-D3 method (Table S1 in ESI), we employed the B3PW91 functional without the D3 dispersion correction in this work. 36


Page 36 of 55

Page 37 of 55 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


(63) This means that two isomers MS2A and MS2B would not arise from the difference in interaction between SO2 and ligand L trans to SO2 . (64) In gas and aqueous phases, we checked the connectivity of GS, MS1, MS2, and TSs by the IRC calculation. In crystal phase, we did not perform the IRC by the large QM/MM calculation (Scheme 3b) in which all QM0, QM1, and QM2 were taken as the QM system, because the computational cost is too large. We instead determined the initial geometry of QM0 part by the small QM/MM calculation in which QM0 was taken as the QM system and the remaining part of crystal was taken as the MM system. The connectivity of these geometries corresponding to GS, MS1, MS2, and TSs was checked using the IRC calculation at the small QM/MM level. It is likely that the connectivity of GS, MS1, MS2, and TSs reasonably remains in the large QM/MM calculation, because the geometry of QM0 does not differ so much between the small and large QM/MM calculations. (65) Su, P.; Li, H. Energy decomposition analysis of covalent bonds and intermolecular interactions. J. Chem. Phys. 2009, 131, 014102. (66) Glendening, E. D.; Badenhoop, J. K.; Reed, A. E.; Carpenter, J. E.; Bohmann, J. A.; Morales, C. M.; Weinhold, F. NBO 5.0. Theoretical Chemistry Institute, University of Wisconsin, Madison, 2001. (67) Consistent with these changes, the Cl negative charge decreases in 1Cl and the positive charges of W and Py increase in 1W and 1Py because charge transfer (CT) from these ligands to the Ru becomes stronger in TSGM than in GS. (68) The weakening of the σ -donation from NH3 to Ru increases the Ru positive charge moderately in the cases of L = W and Py and somewhat in the case of L = Cl− despite of the weakening of the π -back donation from Ru to SO2 . (69) One possible reason is that the out-of-plane non-bonding orbital (nout ) of SO2 forms coordi-

37



nation bond with the Ru in TSMM but the in-plane non-bonding orbital (nin ) forms it in MS1 because the O1 atom of SO2 coordinates with the Ru in MS1; see Scheme 5 (right-hand side). (70) The isomerization from MS2 to MS1 is the most endothermic in 1Cl than in 1W and 1Py , suggesting that TSMM is more product-like in 1Cl than in the others. The Ru-O1 distance increases at TSMM then decreases at MS1. The Ru-O1 of TSMM is the closest to that of MS1 in 1Cl , which shows that TSMM of 1Cl is the most product-like. (71) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999-3093. (72) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal solvation model based on solute electron density and a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J. Phys. Chem. B 2009 113, 6378-6396. (73) This relation, which has been often observed in the calculations based on simulations using explicit solvent molecules, 74,75 is used in the implicit solvation models such as PCM method. (74) Jean-Charles, A.; Nicholls, A.; Sharp, K.; Honig, B.; Tempczyk, A.; Hendrickson, T. F.; Still, W. C. Electrostatic contributions to solvation energies: comparison of free energy perturbation and continuum calculations. J. Am. Chem. Soc. 1991, 113, 1454-1455. (75) Nina, M.; Beglov, D.; Roux, B. Atomic radii for continuum electrostatics calculations based on molecular dynamics free energy simulations. J. Phys. Chem. B 1997, 101, 5239-5248. (76) Two important factors must be considered to explain the difference in CT from NH3 to Ru between GS and MS1. One is the CT from the lσ HOMO orbital of SO2 to Ru dσ , which raises the Ru dz2 orbital energy in GS because of the orbital mixing. The other is the approach of the negatively charged O1 atom of SO2 to the Ru in MS1, which raises the Ru dz2 and dx2 −y2 orbital energy in MS1 due to the electrostatic repulsion with the negatively charged O1 atom. Though the former factor suppresses the CT from the NH3 to Ru in GS but the latter factor

38


Page 38 of 55

Page 39 of 55 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


suppresses the CT in MS1. Tables 1, S6, and S7 in ESI indicate that the CT from the NH3 to Ru is larger in GS than in MS1, which indicates that the latter factor is stronger than the former one in GS and MS1. The other possible reason for the weaker CT in MS1 than in GS is that the approach of the negatively charged O2 atom of SO2 to H atom of NH3 enhances the polarization of the NH3 ligand to lower the lone-pair orbital energy of NH3 , and suppress the CT from the NH3 to Ru in MS1. Actually, the NH3 ligand which is the closest to the O1 atom exhibits the smallest CT in MS1. However, the remaining three NH3 ligands also exhibit the smaller CT in MS1 than in GS, indicating that this is not a major reason to suppress the CT in MS1 relative to GS. (77) Aono, S.; Mori, T.; Sakaki, S. 3D-RISM-MP2 approach to hydration structure of Pt(II) and Pd(II) complexes: Unusual H-ahead mode vs usual O-ahead one. J. Chem. Theory Comput. 2016, 12, 1189-1206. (78) In [PtII (NH3 )4 ]2+ (H2 O), the solvation change (iii) does not occur because the H2 site of W does not form the hydrogen bond with NH3 ligand in the H-interaction form unlike the O2 atom of SO2 in MS1. (79) In the MS2-crystal, the SO2 orientation is the same between target Ru complex and the nearest neighboring Ru complex but in the MS2′ -crystal it differs between them (see Scheme S1 in ESI). (80) The calculation of QM0+QM1+QM2 in gas phase was failed in determining correctly the electronic structure of the target Ru complex QM0, because QM0 was significantly influenced by the counterions QM2. This result clearly indicates the important role of MMcrystals for representing correctly the electronic structure of QM0. Here, the energy component of QM0+QM1+QM2 in the presence of all the MM molecules was taken as the reference in Table 5.

39



Page 40 of 55

Table 1: Change in NBO charge on path II in gas phase.a moiety

GSB

TSGM B

SO2 0.09 -0.10 (-0.19)c Cl -0.49 -0.38 (+0.11) 4Amb 1.27 1.21 (-0.06) Ru 0.13 0.27 (+0.14) 4dσ populations 0.79 0.77 (-0.02) dx2 −y2 dz2 1.00 0.89 (-0.11) SO2 0.28 Wb 0.17 4Amb 1.28 Ru 0.27 4dσ populations dx2 −y2 0.79 dz2 0.85

0.09 (-0.19) 0.21 (+0.04) 1.23 (-0.05) 0.47 (+0.20)

SO2 0.24 Pyb 0.28 4Amb 1.27 Ru 0.21 4dσ populations dx2 −y2 0.78 dz2 0.91

0.77 (-0.02) 0.64 (-0.21)

MS2B 1Cl -0.22 (-0.31)c -0.41 (+0.08) 1.28 (+0.01) 0.35 (+0.22) 0.80 (+0.01) 0.92 (-0.08) 1W -0.04 (-0.32) 0.22 (+0.05) 1.29 (+0.01) 0.53 (+0.26)

TSMM B

MS1B

-0.12 (-0.21)c -0.44 (+0.05) 1.15 (-0.12) 0.41 (+0.28)

-0.17 (-0.26)c -0.46 (+0.03) 1.17 (-0.10) 0.46 (+0.33)

0.74 (-0.05) 0.78 (-0.22)

0.75 (-0.04) 0.79 (-0.21)

0.08 (-0.20) 0.20 (+0.03) 1.17 (-0.11) 0.55 (+0.28)

0.06 (-0.22) 0.18 (+0.01) 1.18 (-0.10) 0.58 (+0.31)

0.73 (-0.06) 0.56 (-0.29)

0.73 (-0.06) 0.58 (-0.27)

0.06 (-0.18) 0.37 (+0.09) 1.21 (-0.06) 0.36 (+0.15)

0.80 (+0.01) 0.71 (-0.14) 1Py -0.06 (-0.30) 0.34 (+0.06) 1.28 (+0.01) 0.44 (+0.23)

0.07 (-0.17) 0.32 (+0.04) 1.16 (-0.11) 0.45 (+0.24)

0.02 (-0.22) 0.31 (+0.03) 1.18 (-0.09) 0.49 (+0.28)

0.75 (-0.03) 0.78 (-0.13)

0.79 (+0.01) 0.81 (-0.10)

0.73 (-0.05) 0.69 (-0.22)

0.73 (-0.05) 0.69 (-0.22)

at the B3PW91 level. Unit is in e. b Am, W, and Py represent ammine, water, and pyridine ligands, respectively. c In parentheses are the changes from that of GSB. a Calculated

40


Page 41 of 55 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


Table 2: Free energy change ∆A in aqueous phase and difference from potential energy change ∆Egas in gas phase.a

GS

TSGM

MS2

∆A ∆A - ∆Egas

-0.1 (-0.4)b

32.0 (+7.0)b

11.1 (+5.7)b

∆A ∆A - ∆Egas

-0.3 (-1.7)

32.4 (+4.0)

11.7 (+2.6)

∆A ∆A - ∆Egas

0.1 (-1.5)

29.0 (+6.5)

9.4 (+3.2)

∆A ∆A - ∆Egas

0.0

31.3 (+5.6b )

11.3 (+4.5)

∆A ∆A - ∆Egas

0.0

32.1 (+3.2)

11.1 (+1.4)

∆A ∆A - ∆Egas

0.0

28.6 (+3.5)

10.0 (+3.7)

TSMM (A) path I 1Cl 37.2 (+11.2)b 1W 31.2 (+17.1) 1Py 25.6 (+15.1) (B) path II 1Cl 37.9 (+11.4) 1W 31.8 (+17.6) 1Py 24.7 (+14.2)

MS1

Barrier TSGM TSMM (MS2→GS) (MS1→MS2)

25.7 (+10.8)b

20.9 (+1.3)b

11.5 (+0.4)b

21.0 (+16.7)

20.7 (+1.4)

10.2 (+0.1)

16.3 (+13.6)

19.6 (+3.3)

9.3 (+1.5)

27.2 (+11.9)

20.0 (+1.1)

10.7 (-0.5)

20.7 (+16.5)

21.0 (+1.8)

11.1 (+1.1)

16.6 (+14.4)

18.6 (+1.8)

8.1 (-0.2)

at the CCSD(T) level. Unit is in kcal mol−1 . GSB on path II is taken as the energy zero. b In parentheses are the difference from the energy in gas phase. a Calculated

41



Page 42 of 55

Table 3: Free energy components on path II in aqueous phase.a

∆A ∆Esolute ∆∆µ ∆Eint /2 Component [SO2 ] [RuII (NH3 )4 L] ∆A ∆Esolute ∆∆µ ∆Eint /2 Component [SO2 ] [RuII (NH3 )4 L] ∆A ∆Esolute ∆∆µ ∆Eint /2 Component [SO2 ] [RuII (NH3 )4 L]

GSB

TSGM B

MS2B

0.0 (0.0b ) (0.0) (0.0) (0.0)

31.3 (30.3) (23.6) (6.7) (6.9)

11.3 (12.3) (6.7) (5.6) (6.0)

(0.0) (0.0)

(5.0) (1.9)

(8.6) (-2.6)

0.0 (0.0) (0.0) (0.0) (0.0)

32.1 (32.2) (29.6) (2.6) (2.4)

11.1 (11.1) (8.9) (2.2) (2.3)

(0.0) (0.0)

(6.5) (-4.1)

(14.6) (-12.3)

0.0 (0.0) (0.0) (0.0) (0.0)

28.6 (28.4) (22.8) (5.6) (5.2)

(0.0) (0.0)

(9.8) (-4.6)

TSMM B MS1B 1Cl 37.9 (34.0) (22.2) (11.8) (10.2) (5.8) (4.4) 1W 31.8 (31.3) (11.2) (20.1) (18.4)

Barrier TSGM B TSMM B (MS2B→GSB) (MS1B→MS2B)

27.2 (23.2) (13.7) (9.5) (7.6)

20.0 (18.0) (16.9) (1.1) (0.9)

10.7 (10.8) (8.5) (2.3) (2.6)

(6.6) (1.0)

(-3.6) (4.5)

(-0.8) (3.4)

20.7 (18.3) (-0.9) (19.2) (18.3)

21.0 (21.1) (20.7) (0.4) (0.1)

11.1 (13.0) (12.1) (0.9) (0.1)

(5.1) (13.2)

(-8.1) (8.2)

(-2.1) (2.2)

10.0 (11.5) (6.4) (5.1) (5.4)

(3.0) (15.4) 1Py 24.7 (25.7) (7.4) (18.3) (17.7)

16.6 (14.5) (-2.5) (17.0) (16.7)

18.6 (16.9) (16.4) (0.5) (-0.2)

8.1 (11.2) (9.9) (1.3) (1.0)

(17.4) (-12.0)

(5.6) (12.1)

(6.3) (10.4)

(-7.6) (7.4)

(-0.7) (1.7)

at the CCSD(T) level. Unit is in kcal mol−1 . GSB on path II is taken as the energy zero. b In parentheses are calculated at the B3PW91 level. a Calculated

42


Page 43 of 55 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


Table 4: Potential energy change ∆Egas in gas phase, free energy change ∆A in aqueous phase, and potential energy change Ecry in crystal phase by thermal isomerization of 1Py .a

∆Egas ∆A ∆Ecry

path I path II path I path II MS1-crystal 100 Kc 120 K 140 K MS2-crystal 100 Kc 120 K 140 K

GS

TSGM

1.6 0.0b 0.1 0.0b

Barrier TSGM TSMM (MS2→GS) (MS1→MS2) 16.3 7.8 16.8 8.3 19.6 9.3 18.6 8.1

MS2 TSMM

MS1

22.5 23.1 29.0 28.6

6.2 6.3 9.4 10.0

10.5 10.5 25.6 24.7

2.7 2.2 16.3 16.6

(0.0)b (0.0) (0.0)

31.7 31.7 32.2

15.2 15.2 15.3

38.0 38.1 39.9

27.1 28.0 31.3

16.5 16.5 16.9

10.9d 10.1 8.6

0.0b 0.0 0.0

31.7 31.8 32.4

14.6 14.7 15.3

39.9 40.5 42.0

29.7 30.5 32.7

17.1d 17.1 17.1

10.1 10.0 9.3

at the CCSD(T) level. Unit is in kcal mol−1 . Geometry was determined in each of gas, solution, and crystal phases. b Reference values. c Lattice vector at 100, 120, and 140 K, respectively. d MS1→MS2 isomerization occurs in the MS1-crystal and the MS2→GS isomerization occurs in the MS2-crystal. a Calculated

43



Page 44 of 55

Table 5: Influences of surrounding QM and MM species on target QM molecule QM0 in 1Py .a

QM

MM

GS

QM0b 0+1b 0+2b 0+1+2c 0+1+2c 0+1+2c

none none none none MM1 to 8 all MM

0.0 0.0 0.0 0.0 0.0 0.0

QM0b 0+1b 0+2b 0+1+2c 0+1+2c 0+1+2c

none none none none MM1 to 8 all MM

0.0 0.0 0.0 0.0 0.0 0.0

TSGM

MS2

TSMM

MS1

MS1-crystal at the lattice vector at 100 K 24.3 4.7 8.3 1.8 27.2 (+2.9)d 7.7 (+3.0)d 18.4 (+10.1)d 14.5 (+12.7)d 30.1 (+5.8) 16.6 (+11.9) 22.9 (+14.6) 15.8 (+14.0) 36.0 (11.7) 22.3 (17.6) 41.6 (+33.3) 32.5 (+30.7) 34.1 [-1.9] 20.0 [-2.3]e 41.8 [+0.2]e 30.5 [-2.8]e 31.7 [-2.4] 15.2 [-4.8] 38.0 [-3.8] 27.1 [-3.4] MS1-crystal at the lattice vector at 140 K 24.7 5.8 9.8 3.3 26.8 (+2.1)d 7.9 (+2.1)d 19.2 (+9.4)d 19.6 (+16.3)d 28.2 (+3.5) 14.8 (+9.0) 22.7 (+12.9) 17.9 (+14.6) 33.8 (+9.1) 20.5 (+14.7) 42.4 (+32.6) 41.4 (+38.1) 31.2 [-2.6] 18.4 [-2.1]e 39.5 [-2.9]e 36.6 [-4.8]e 32.2 [+1.0] 15.3 [-3.1] 39.9 [+0.4] 31.3 [-5.3]

a Calculated

Barrier TSGM TSMM (MS2→GS) (MS1→MS2) 19.6 19.5 (-0.1)d 14.5 (-5.1) 13.7 (-5.9) 14.1 [+0.4]e 16.5 [+2.4]

7.0 3.9 (-3.1)d 7.1 (+0.1) 9.1 (+2.1) 11.3 [+2.2]e 10.9 [-0.4]

18.9 18.9 (+0.0)d 13.4 (-5.5) 13.3 (-5.6) 12.8 [-0.5]e 16.9 [+4.1]

6.5 -0.4 (-6.9)d 4.8 (-1.7) 1.0 (-5.5) 2.9 [+1.9]e 8.6 [+5.7]

at the geometry determined in the MS1-crystal. Self energy of target molecule was evaluated at the CCSD(T) level, while interaction energies with other QM molecules were evaluated at the B3PW91 level. Unit is in kcal mol−1 . QM0 represents the target Ru complex 1Py . QM1 represents two neighboring Ru complexes. QM2 represents one neighboring water molecule and five neighboring Cl− counterions. MM1 to MM8 represent the nearest neighboring eight MM molecules; see Scheme S1 for their positions. b QM calculation was performed in the absence of the MM molecules in MS1-crystal. c QM calculation was performed in the presence of all MM molecules in MS1-crystal. d In parentheses are the difference from the first column evaluated for QM0. e In brackets are the difference from one column above.

44


Page 45 of 55 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


Isomerization of Sulfur-Dioxide Comlex O1 NH3

O2 S

NH3 hν NH3

Ru NH3

O1 S

O2 NH3 hν NH3

Ru NH3

NH3 L

GS ground state η1 -structure

∆

NH3

NH3 L

∆

MS2 meta-stable η2 -structure

S O2 O NH3 Ru NH3 L 1

MS1 meta-stable η1 -structure

z y

x

1a: L = Cl − 1b: L = H2O 1c: L = pyridine

Scheme 1: Isomerization of ruthenium(II) sulfur dioxide complex in crystal.

45



unoccupied σ*

Page 46 of 55

* π out

* π in

LUMO

antibonding 6.22 eV (2.93 eV)

Occupied HOMO lσ

5.72 eV (2.33 eV)

n in

0.87 eV (-3.85 eV)

n out

nonbonding -13.50 eV (-9.37 eV)

-14.90 eV (-10.25 eV)

π in

σ

-14.25 eV (-10.62 eV)

π out

bonding -19.10 eV (-13.82 eV)

-18.81 eV (-13.70 eV)

In-plane π

σ

-18.02 eV (-13.77 eV)

Out-of-plane π

Scheme 2: In-plane and out-of-plane π -orbitals of SO2 .a a Orbital energy was evaluated with Hartree-Fock method using the same basis sets described in computational section. In parentheses are evaluated at the R-DFT/B3PW91 level.

46


Page 47 of 55 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) Construction of periodic MM crystal (GS-, MS1-, and MS2-crystal). MM crystal is assumed. Φ I is solved by Fock matrix derived from Eq. (2).

SCF (charge)

QM ESP charge is evaluated. SCF(charge) is completed. I

Energy gradient of Ecrystal is evaluated (Ref. [32]). New geometry of QM subunit

SCF (geometry)

SCF(geometry) is completed.

Periodic MM crystal is determined. (b) Optimization of TS and refinement of GS, MS1, and MS2 in MM crystal.

Φ is calculated for target QM under fixed MM crystal (Eq. (7)). Energy gradient is evaluated (Eq. (9)). New geometry of target QM molecule.

Converged.

not converged.

Scheme 3: Flowchart of the present QM/MM procedure to investigate reaction of target QM molecule in crystal.

47



Page 48 of 55

(a) Construction of periodic MM crystal.

Target Ru complex H2O Cl

Cl

SO2 Ru

Py

NH3

NH3

SO2

Cl

the nearest neighboring pyridine

Ru Py H2O Cl

z

one neighboring Ru complex

y

(b) Optimization of TS and refinement of GS, MS1, and MS2

Target Ru complex Ru

QM1

QM2 five neighboring Cl neighboring H 2 O

SO2 NH3

Cl Py

Ru SO2

Py

the nearest neighboring pyridine

SO2

QM0

NH3

H2O z

Ru y

Py

two neighboring Ru complexes (QM1)

QM1

Scheme 4: QM region in the QM/MM calculation for (a) construction of periodic MM crystal and (b) optimization of TS and refinement of GS, MS1, and MS2.

48


Page 49 of 55 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) path I

S

S

O1

O2

1

O Ru O2

NH3

MS2A

GSA

2 S O

O1 MS1A

2

O is close to one NH 3 and S exists between two NH3 . SO 2 is eclipsed to Py plane.

(b) path II O1

O1

S

Ru NH3

2

O

GSB

O2 MS2B

O1 S S 2

O MS1B

2

O exists between two NH3 and S is close to one NH 3 . SO 2 is staggered to Py plane.

Scheme 5: Difference in geometry between paths I and II.

49



(A) Solvation of [RuII (NH3 )4 (SO2 )L]n+ In GS

In MS1

change-(i) solvation space around O1 is small but g(H) increases g(H) is large. around 4Am. region b region a change-(ii) region c solvation space region c around S S O2 O2 O1 increases. region a S δ+ δ- O1 change-(iii) region b solvation space Am Am Am Am around O2 2+ 2+ decreases. δ-: CT M M solvation space around S is small but g(O) increases around 4Am.

Am

region d

g(O) is large.

Am

Am

L

Am

L

change-(iv) g(O) decreases around 4Am and L.

region d

(B) Comparison between [RuII (NH3 )4 (SO2 )L]n+ and [PtII (NH3 )4 ]2+ -W O1

O2

δ+ S

Am

Am Ru

Am

L

Am

S δ+(S) stabilizes d x2−y 2.

δ-

Am

L

δ − (O1) destabilizes d x2−y 2.

δ-: CT CT(4Am Ru) decreases.

Am

MS1

O Am

Pt Am

Am Ru

H2

Am

O1

Am

CT(4Am Ru) is large.

GS

δ+H1

O2

Am

H2

H1 δ+(H1) stabilizes d x2−y 2. CT(4Am Pt) is large.

H-interaction form

O

δ-

Am

Am Pt

Am

δ − (O) destabilizes d x2−y 2.

δ-: CT CT(4Am Pt) decreases. Am

O-interaction form

Scheme 6: (A) Difference in solvation between GS and MS1 and (B) comparison of charge re-distribution between [RuII (NH3 )4 (SO2 )L]n+ and [PtII (NH3 )4 ]2+ -W.

50


Page 50 of 55

Page 51 of 55 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


top view QM1 QM2

O1 Cl

Py

QM1 QM2

H Py

O1

2

Ru

S O2 O is close to H Py .

Ru

QM0

QM1 Py

Cl

Cl S

H Py

O1

O 2 leases O2 H . Py

S

Py

QM2

Ru

H Py O2 S is close to H Py .

QM0

QM0

side view QM2Cl S is close

QM2 Cl

to Cl.

QM0

S

Py

QM1 (A) GS

QM2 Cl O 1 is close

O 2 is close O2 to Cl.

QM0

Py

to Cl.

QM0

QM1 (B) MS2

O1

Py

QM1

(C) MS1

Scheme 7: Characteristic difference in intra- and inter-molecular interactions among GS, MS2, and MS1 structures.

51



Page 52 of 55

(a) path I L = pryridine

3.15

O1

S

2 O2 S O O1 S O2 1 2.18 2.86 2.50 2.17 O 2.52

Ru

Ru

1 2.31 O

Ru

S S O2 1 2.11 O 3.16

Ru

O2 3.35

Ru

L

L

L

L

GSA

TSGM A

MS2A

TSMMA

MS1A

δr = -0.98

δr = -0.36

δr = 0.35

δr = 0.85

δr = 1.24

L

(b) path II L = pryridine O1 O2 S 2.17 3.15 Ru

S 3.38 S O2 O1 S O2 2 1 O 2.28 1 O 3.07 S 2.47 1 O2 2.85 2.50 2.18 O 2.12 O Ru

Ru

Ru

Ru

L

L

L

L

GSB

TSGM B

MS2B

TSMMB

MS1B

δr = -0.97

δr = -0.35

δr = 0.31

δr = 0.79

δr = 1.26

L

Figure 1: Changes in geometry of ruthenium(II) sulfur dioxide complex with L = pyridine in gas phase; see Figure S1 in ESI for comparison with those in aqueous phase and MS1-crystal at 140 K.

52


Page 53 of 55

(A) In gas phase (a) path I

Egas [kcal/mol]

30

Eǂ

28.4 (19.6) 25.0 (19.3) 26.0 22.5 (16.3)

25 20

1W 1Cl (b) path II

Eǂ

30

(11.1) (10.1) (7.8)

15

9.1 6.2 5.4

10 5

1.6 1.4 0.3

14.1 10.5

20

4.0 2.7

5

0.1 -0.1 -0.3

29.0

31.2 25.6

Aǂ

11.7

(20.9) 11.0 (20.7) 9.4 (19.6)

Aǂ

0.0 0.0 0.0

14.2 10.5

Eǂ

(11.2) (10.0) (8.3) 15.3 4.2 2.2

GSB TSGM B MS2B TSMM B MS1B

(B) In aqueous phase (c) path I 37.2

9.7 6.8 6.3

10

0

32.4 32.0

Eǂ

28.9 (18.9) 25.7 (19.2) 26.5 23.1 (16.8)

25

GSA TSGM A MS2A TSMM A MS1A

40 35 30 25 20 15 10 5 0 -5

1Py

15

14.9

0

A [kcal/mol]

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


40 35 30 25 20 15 10 5 0 -5

25.7 21.0 16.3

(11.5) (10.2) (9.3)

GSA TSGM A MS2A TSMM A MS1A

(d) path II 32.1 31.3

37.9 31.3

28.6

24.7

0.0 0.0 0.0

Aǂ

11.3 11.1

(20.0) 10.0 (21.0) (18.6)

Aǂ

27.2 20.7 16.6

(10.7) (11.1) (8.1)

GSB TSGM B MS2B TSMM B MS1B

Figure 2: Changes in potential energy by thermal isomerization at the CCSD(T) level in gas phase ((a) on path I and (b) on path II) and that in free energy of thermal isomerization at the CCSD(T) level in aqueous phase ((c) on path I and (d) on path II).

53



Lattice vector: 100 K

120 K

Ecrystal [kcal/mol]

45 40 35 30 25 20 15 10 5 0

43.5 41.8 32.2 41.6 29.0 18.4 28.8 17.5 Eǂ Eǂ 17.4 (11.3) (16.7) 35.1 33.7 33.5

0.0 (12.8) 0.0 (16.3) 0.0 (16.0) (12.8) GS TSGM MS2 TSMM MS1

(c) MS2-crystal 32.4 31.8 31.7

140 K

(b) MS1-crystal

(a) GS-crystal 45 40 35 30 25 20 15 10 5 0

Ecrystal [kcal/mol]

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 54 of 55

42.0 32.7 40.5 39.8 30.5 29.7

E ǂ 15.3 Eǂ (9.3) 0.0 (17.1) 14.7 14.6 (17.1) (10.0) 0.0 0.0 (17.1) (10.1) GS TSGM MS2 TSMM MS1

45 40 35 30 25 20 15 10 5 0

45 40 35 30 25 20 15 10 5 0

32.2

39.9 31.3 38.1 31.7 38.0 31.7 28.0 27.1 15.3 Eǂ Eǂ 15.2 0.0 (16.9) 15.2 (8.6) (16.5) (10.1) 0.0 0.0 (16.5) (10.9) GS TSGM MS2 TSMM MS1

(d) MS2'-crystal 32.6 31.9 31.8 E ǂ 15.2 0.0 (17.4) 14.7 0.0 (17.2) 14.6 0.0 (17.2) GS TSGM MS2

Figure 3: Changes in potential energy by thermal isomerization of [Ru(NH3 )4 (SO2 )Py](2Cl·H2 O) at the CCSD(T) level in crystal.

54


41.7 32.8 40.6 39.6 30.7 29.6 Eǂ

(8.9) (9.9) (10.0) TSMM MS1

Page 55 of 55

Graphical TOC Entry ∆

∆

S O2 O1 Ru

∆E or ∆A

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


O1 S

crystal

O2

Ru O1

solution

O2 S

Ru

GS

gas TS GM

MS2

TS MM

MS1

55


MM Approach to Isomerization of Ruthenium(II) Sulfur Dioxide

Recommend Documents