Coupled Cluster Theory Combined with Reference Interaction Site

Mar 29, 2018 - Coupled Cluster Theory Combined with Reference Interaction Site Model Self-Consistent Field Explicitly Including Spatial Electron Densi...
0 downloads 3 Views 276KB Size
Subscriber access provided by University of Florida | Smathers Libraries

Condensed Matter, Interfaces, and Materials

Coupled cluster theory combined with reference interaction site model selfconsistent field explicitly including spatial electron density distribution Daisuke Yokogawa J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00168 • Publication Date (Web): 29 Mar 2018 Downloaded from http://pubs.acs.org on March 30, 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 19 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

Journal of Chemical Theory and Computation

Coupled cluster theory combined with reference interaction site model self-consistent field explicitly including spatial electron density distribution D. Yokogawa∗,†,‡ †Department of Chemistry, Graduate School of Science, Nagoya University, Chikusa, Nagoya 464-8602, Japan ‡Institute of Transformative Bio-Molecules (WPI-ITbM), Nagoya University, Chikusa, Nagoya 464-8602, Japan E-mail: [email protected] Phone: +81 52 789 2851 Abstract Calculating the geometry and energy of a molecule in a solution is one of the most important tasks in chemistry. However, performing an accurate calculation in a solution is still a difficult task because the electronic structure and solute-solvent interactions are required to be accurately evaluated with an efficient computational cost. To overcome this difficulty, we proposed the coupled-cluster with single and double excitations and perturbative triple excitations combined with the reference interaction site model (RISM) by employing our fitting approach. Our method correctly reproduced the relative stabilities of 1,2,3-triazole, isonicotinic acid, cytosine, and 6-chloro-2-pyridone in the aqueous phase, whereas the dielectric continuum model provided incorrect results

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

for isonicotinic acid and 6-chloro-2-pyridone. Our method provided accurate results because the RISM captured the local solvation structure, such as hydrogen bonds.

1

Introduction

Calculating the geometry and energy of a molecule in a solution is one of the most important tasks in chemistry. An accurate prediction of the geometry and energy can accelerate the understanding of chemical reaction mechanisms in a solution for several fields, including organic chemistry, photophysical chemistry, and biochemistry. To date, many solvation theories have been proposed. 1–3 However, performing an accurate calculation in a solution is still a difficult task because the electronic structure and solutesolvent interactions are required to be accurately evaluated with an efficient computational cost. In the gas phase, coupled-cluster with single and double excitations and perturbative triple excitations (CCSD(T)) is well known as the gold standard for quantum chemistry and yields reliable results. 4–6 However, to the best of our knowledge, there is no such universally recognized gold standard for solvation theory. The polarizable continuum model (PCM) is a popular solvation theory. The equations employed in PCM are simple, and the computational cost is low. As such, PCM has been combined with coupled-cluster (CC) methods. 7,8 Although this combination is commonly employed, there is an inherent weakness in this strategy: the computed local solvation interactions, such as hydrogen bonding, are underestimated. The reference interaction site model (RISM) is another solvation theory, and it has the potential to overcome such PCM defects. 9,10 RISM has been coupled with quantum mechanics (QM), which is referred to as RISM self-consistent field (RISM-SCF), and applied to certain systems. 11–13 Although RISM-SCF can accurately evaluate hydrogen bonding, it has another severe weakness: the calculations are unstable in polar solvents. Even for a simple system, for example, a system involving methanol, we cannot obtain the converged results with the original RISM-SCF

2

ACS Paragon Plus Environment

Page 2 of 19

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

Journal of Chemical Theory and Computation

theory. 14 To overcome the limitation of PCM and the original RISM-SCF theory, we developed a new method based on the RISM and our new fitting approach, which is referred to as RISM-SCF-SEDD. By including spatial electron density distribution, the instability problem associated with RISM-SCF was resolved to a considerable extent. 14–16 Herein, we coupled the RISM-SCF-SEDD method with CCSD(T) and demonstrated the reliability of this new method. The context of this paper is as follows. In Section 2, we derive the CC equations coupled with RISM. The computational details are summarized in Section 3. To discuss the accuracy of the method, four tautomeric molecules were selected for performing calculations in a solution. The basis set dependency and the relative stabilities are summarized in Section 4.

2 2.1

Method Derivation of CC equations coupled with RISM

In RISM, we compute correlation functions between particles. The total correlation function hαs (r) and direct correlation function cαs (r) between the solute α site and the solvent s site are computed with the following equations: ∑

[ωαγ ∗ cγt ∗ χts ] (r),

(1)

hαs (r) = exp [−βuαs (r) + hαs (r) − cαs (r)] − 1,

(2)

hαs (r) =

γt

where ωαγ is the intramolecular correlation function of the solute molecule, χts is the pure solvent site density pair correlation function, uαs is the site-site pair potential, and the asterisk denotes a convolution integral. The radial distribution function (RDF) between the solute α site and the solvent s site is obtained with hαs (r) + 1. Equation 2 is called as

3

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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 4 of 19

hypernetted-chain (HNC) closure. The free energy derived from the HNC closure is [ ] ∫ 1∑ 1 2 1 ∆µ = ρs dr hαs (r) − cαs (r) − hαs (r)cαs (r) . β αs 2 2

(3)

hαs and cαs are computed with eqs 1 and 2 under the site-site pair potential uαs that is defined with,

uαs (r) =

unES αs (r)



∑ i∈α

∫ di

fi (r′ ) ′ dr , |r − r′ |

(4)

where unES αs is non-electrostatic site-site pair potential and {fi } is the auxiliary basis set (ABS). {di } is the fitting coefficient and determined by

d=



(Ξ + Γ)−1 R′pq γpq ,

(5)

pq ′

where p and q labels are employed for the general MOs. Ξ, Γ, and Rpq are matrices that have been defined in previous papers, 15,17 and γpq is the reduced one electron density matrix. CC In the case of CC, one electron density matrix γpq is given by, 18

⟩ ⟨ CC γpq = Φ0 (1 + Λ)e−T {a†p qq }eT Φ0 ,

(6)

where Φ0 is the reference function computed in a solution, T is the excitation operator, {λ} CC are de-excitation amplitudes, and Λ is the de-excitation operator. If we employ the γpq

in eq 5, hαs and cαs depend on Φ0 , T , and {λ}, and ∆µ becomes functional of these three functions, i.e. ∆µ[Φ0 , T, Λ]. As discussed in the coupling between CC and PCM, decoupling the CC equations from the Λ equations is the key point for proposing an efficient method. 8 Herein, to derive the

4

ACS Paragon Plus Environment

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

Journal of Chemical Theory and Computation

equations, we defined the following Lagrangian: ∑ ∂∆µ 0⟩ −T T 0 λµ , L = Φ (1 + Λ)e H0 e Φ + ∆µ[Φ , T, 0] + ∂λ µ 0 µ ⟨

0

(7)

where H0 is the molecular Hamiltonian. The derivative of the solvation free energy, ∆µ, with respect to the λ amplitude at λ = 0 is expressed as follows: 19 ∑[ ⟩ ] ⟨ −T † ∂∆µ t −1 ′ e {ap qq }eT Φ0 . = V (Ξ + (m − 1)Γ) R µ pq ∂λµ 0 pq

(8)

The CC equation can be derived by minimizing L with respect to λ: ⟩ ∂∆µ ⟨ −T ∂L T 0 = µ e H0 e Φ + ∂λµ ∂λµ 0 ⟩ ⟨ = µ e−T H0solv eT Φ0 (9)

= 0,

where H0solv is the solvated Hamiltonian, which is expressed as follows: H0solv ≡ H0 +

∑[ ] Vt (Ξ + (m − 1)Γ)−1 R′pq {a†p aq }.

(10)

pq

By determining T amplitudes using eq 9, we can obtain the CC free energy as follows: ⟩ ⟨ T (Λ) GCC = Φ0 e−T H0 eT Φ0 + ∆µ[Φ0 , T, 0],

(11)

where superscript T (Λ) indicates that we considered T amplitude and λ amplitude up to T (Λ)

the first order in the derivation of the free energy. By adding (T) correction to GCC , we defined the CCSD(T) free energy in a solution as follows: 5,20 T (Λ)

T (Λ)

[4]

[5]

GCC(T) = GCC + ET + EST .

5

ACS Paragon Plus Environment

(12)

Journal of Chemical Theory and Computation 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

2.2

Page 6 of 19

Comparison with the previous approach

In previous studies, 21–24 we defined CCSD free energy in a solution as follows: ⟩ ⟨ G0CC = Φ0 e−T H0 eT Φ0 + ∆µ[Φ0 , 0, 0].

(13)

The superscript 0 in G0CC indicates that the solvation structure is determined with the HF wave function and does not depend on T and Λ. The CCSD(T) free energy was computed using the following equation: [4]

[5]

G0CC(T) = G0CC + ET + EST .

(14)

The difference between the previous research and this study resides in the inclusion or exclusion of the T amplitude in ∆µ. Because ∆µ in the previous study did not depend on the T amplitude, the correction in the CC calculations did not affect the RISM results and the solvation structures were calculated using the HF wave function. In the proposed method, the solvation effect is determined with the HF wave function and T amplitude and the solvation structure is updated during the CCSD iteration. As discussed later, the additional computational cost of including the T amplitude in ∆µ is not large because the computational cost of RISM is much smaller than that of CCSD(T).

2.3

Linear dependency in the auxiliary basis set (ABS)

In the RISM-SCF-SEDD, spherical Gaussian functions were employed for ABS preparation and the exponents were automatically determined based on the exponents of the basis set employed in the QM calculation. When we employ a large basis set in the QM calculation, there is a near-linear dependence in the constructed ABS. This occasionally leads to

6

ACS Paragon Plus Environment

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

Journal of Chemical Theory and Computation

numerical problems. To address this issue, we defined a modified ABS, {f˜}, as follows: ˜ ˜ , f = Uf

(15)

where U is defined with 

˜ (1) 0 ··· 0  U   0 U ˜ (2) · · · 0 ˜ = U  . .. ..  .. . .   ˜ (M ) 0 0 ··· U

     ,   

(16)

and  ˜ (I)

U

  =  

 (I) U1,1

(I) U1,2

.. .

.. .

(I)

(I)

···

(I) U1,K−m

.. . (I)

UK,1 UK,2 · · · UK,K−m

  .  

(17)

Here, K is the number of original ABSs on the I site and m is the number of the functions removed from the original ABS. To determine the removed functions, we diagonalized the block matrix of Ξ as follows:

U(I)t Ξ(I) U(I) = ξ (I) .

(18)

The element of the block matrix Ξ(I) is given by ∫∫ (I) Ξij

=

fi (r1 )|r1 − r2 |fj (r2 )dr1 dr2 , (i, j ∈ I).

(19)

When there are small values for ξ (I) , the inverse calculation in eqs 8 and 10 leads to numerical instability. By removing the small eigen values in ξ (I) and reordering U(I) , we constructed ˜ (I) . U

7

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

3

Page 8 of 19

Computational details (a)

(b)

O H

N C

1H

N

2H

Z

(c)

C1

C3

C8A

C8C

(d) Cl

K

EA

EC

Figure 1: Conformations considered for (a) 1,2,3-triazole, (b) isonicotinic acid, (c) cytosine, and (d) 6-chloro-2-pyridone. The label of each conformer is shown under the structure. All calculations were performed in water. Solvation effects were evaluated with RISMSCF-SEDD and PCM coupled with SMD. 25 In RISM calculation, the number density of water was 0.033426 Å−3 and temperature was 300 K. The employed Lennard-Jones parameters are summarized in Supporting information (Table S1). In PCM-CCSD calculation, we employed PTES scheme, where the explicit PCM terms explicitly appear in the iterative solution of the CC equations. 8 We optimized the geometries at the B3LYP level. The basis set for optimization was cc-pVDZ (DZ). 26,27 For CC calculations, we also employed four other basis sets: aug-cc-pVDZ (aDZ), 26–28 cc-pVTZ (TZ), 26,27 and maug-cc-pVTZ (maTZ). 29 Geometry optimization and CC calculation with PCM were performed with the Gaussian 09 and 16 programs, respectively. 30,31 All other calculations were performed with the modified 8

ACS Paragon Plus Environment

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

Journal of Chemical Theory and Computation

GAMESS program, 32 where the RISM-SCF-SEDD code was implemented by us. The conformations considered herein are shown in Figure 1. We focused on the free energy differences among these conformations. Experimentally, the stabilities of cytosine and 6-chloro-2-pyridone were discussed on the basis of pKa and it is considerably difficult to distinguish between the two imine conformers C8A and C8C and the two lactim conformers EA and EC. To compare our results with the experimental data, we computed the effective free energy differences between C1 and the imines ∆Geff C1→C8 and between K and lactim ∆Geff K→E using the following equation: ] 1 [ −β∆GC1→C8A ∆Geff + e−β∆GC1→C8C , C1→C8 = − ln e β ] 1 [ −β∆GK→EA + e−β∆GK→EC . ∆Geff K→E = − ln e β

4

(20) (21)

Results and discussion

The CC basis set dependency in the gas phase has been throughly discussed elsewhere. 33,34 In this study, we focused on the CC basis set dependency with regards to solvation effects. The computed free energy differences are summarized in Table 1. When a small basis set, such as DZ, was used, the difference between the HF and the CCSD(T) was large, whereas the difference was small in the case of maTZ. This is because the large basis set can partially include the dynamic correlation. If a small basis set, such as DZ, is employed, we cannot compute the dynamic correlation accurately, even with CCSD(T). By including an augmented basis set (aDZ), or splitting the valence orbital (TZ), the results become close to those with maTZ. When computer resources are limited and using maTZ becomes prohibitive, aDZ appears to perform better than TZ, though there is difference (less than 1 kcal/mol) between aDZ and maTZ. T (Λ)

In Table 1, we also showed ∆GCC to check the importance of connected triple excitations in these systems. In the previous study about reaction enthalpies, it was mentioned that we

9

ACS Paragon Plus Environment

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

Page 10 of 19

T (Λ)

Table 1: Bais set dependency in the relative free energies, ∆GHF , ∆G0CC(T) , ∆GCC(T) , and T (Λ)

∆GCC . In the cases of cytosine and 6-chloro-2-pyridone, the effective free energy differences computed with eqs 20 and 21 are shown. (Unit: kcal/mol) DZ 1,2,3-Triazole ∆GHF ∆G0CC(T) T (Λ) ∆GCC T (Λ) ∆GCC(T)

aDZ

TZ maTZ

1H → 2H -0.62 0.09 0.08 0.50 -0.21 0.27 0.09 0.51

-0.81 -0.31 -0.60 -0.30

1.18 3.39 3.19 3.21

N→Z -1.03 -2.00 -2.09 -1.43 -1.62 -1.40 -2.17 -1.59

-1.33 -1.33 -1.12 -1.46

2.81 3.25 3.32 3.20

C1 → C3 3.20 3.04 3.70 3.51 3.84 3.71 3.70 3.57

2.88 3.29 3.48 3.33

8.17 6.02 5.40 5.69

C1 → C8 8.44 8.29 7.02 6.63 6.40 6.33 6.93 6.69

7.36 6.04 5.54 6.00

1.18 0.23 0.07 0.26

Isonicotinic acid ∆GHF ∆G0CC(T) T (Λ) ∆GCC T (Λ) ∆GCC(T) Cytocine ∆GHF ∆G0CC(T) T (Λ) ∆GCC T (Λ) ∆GCC(T)

∆GHF ∆G0CC(T) T (Λ) ∆GCC T (Λ) ∆GCC(T)

6-Chloro-2-pyridone ∆GHF ∆G0CC(T) T (Λ) ∆GCC T (Λ) ∆GCC(T)

-0.75 -0.17 0.08 0.02

K→E 0.07 -1.19 0.68 -1.11 0.96 -0.71 0.79 -0.89

10

ACS Paragon Plus Environment

0.63 0.35 0.82 0.53

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

Journal of Chemical Theory and Computation

must take into account the effects of the connected triples. 6 However, at least in our systems, T (Λ)

T (Λ)

the difference between ∆GCC and ∆GCC(T) is less than 0.6 kcal/mol and the trend of the relative stabilities was reproduced with CCSD. T (Λ)

Table 2: Coupled cluster iteration cycles required in the calculations of G0CC and GCC . C1 Z 2H K 22 24 19 22 21 21 17 22

G0CC T (Λ) GCC

In our previous studies, the solvation effect was computed using the HF wave function. Because the CC equation of this study is decoupled from the Λ equations, the CPU time for T (Λ)

∆GCC(T) calculation is only a few percent longer than for ∆G0CC(T) (3.2% and 1.3% longer in the case of C8A and C8C, respectively). Moreover, RISM did not affect the convergence of the CC equations. We summarized the coupled cluster iteration cycles of C1, Z, 2H, T (Λ)

and K required in the calculations of G0CC and GCC in Table 2. From the computational time and the CC iteration cycles, we can say that the computational cost of the present our method does not change so much compared to that of the previous our method. Although the T (Λ)

difference between ∆G0CC(T) and ∆GCC(T) was small (less than 0.4 kcal/mol) in this study, it will become important in the application of RISM to coupled cluster response theories in excited state calculations. 35 Table 3: Relative free energies computed with CCSD(T) using our method (RISM) and PCM. As reference values, the experimental data are also shown. (Unit: kcal/mol) 1H → 2H -3.59 -0.85 -0.30

N→Z 32.65 2.35 -1.46

-0.38,a -2.45b

-2.44 ∼ -2.59c

gas PCM RISM exp.

a

Ref 36 , b Ref 37 ,

c

Ref 38 ,

d

Ref 39 ,

e

C1 → C3 C1 → C8 7.00 1.00 4.32 7.04 3.33 6.00 3.6d

5.5∼6.9d

K→E -3.92 -2.30 0.53 0.46e

Ref 40

The relative free energies computed in the gas and aqueous phases using PCM and RISM 11

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

are summarized in Table 3 alongside experimental values for reference. For 1,2,3-triazole and cytosine, the relative stabilities in solution were similar to those in the gas phase, and PCM could reproduce the experimental data well. On the other hand, the hydration effect changed the most stable conformer of isonicotinic acid and 6-chloro-2-pyridone. PCM could not reproduce the relative stabilities while our method successfully calculated the experimental result. 3 "mACCT_Z.rdf" using 1:27 OH O "mACCT_N.rdf" using 1:27

O

Radial distribution function

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

Page 12 of 19

O

2 N

N H

1

0 0

2

4

6

8

10

r (angstrom)

Figure 2: Radial distribution functions (RDFs) between the solvent hydrogen site and one of the oxygen sites (circled in the inset). The RDFs of the Z and the N forms are shown with solid and dotted lines, respectively. The inability of PCM to calculate isonicotinic acid in solution has already been noted by Nagy et al. 38 Solvent water greatly stabilizes the conformer Z. To compute the solvation energy correctly, Nagy et al. employed a free-energy perturbation method as implemented in Monte Carlo simulations. By considering explicit water molecules, the relative stability of isonicotinic acid was correctly reproduced. Using RISM, although explicit water molecules are not considered, we can calculate hydrogen bonding correctly with an efficient computational cost. In Figure 2, the radial distribution function between the solvent hydrogen site and one of the oxygen sites of isonicotinic acid is shown (circled with dashed line in the inset). The peak around 2 Å in the Z form is higher than that in N, which shows that the hydrogen bonding between solvent water and the Z form is stronger than that between solvent water and the N form. The strong hydrogen bonding greatly stabilizes the Z form 12

ACS Paragon Plus Environment

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

Journal of Chemical Theory and Computation

and RISM can successfully evaluate the hydrogen bondings.

5

Conclusions

We proposed a new method by coupling CCSD(T) with RISM using our fitting approach (SEDD). Although the original RISM-SCF was unstable with large basis sets, RISM-SCFSEDD can perform the coupled cluster calculations with large basis sets, such as cc-pVTZ and maug-cc-pVTZ. Our method correctly reproduced the relative stabilities of 1,2,3-triazole, isonicotinic acid, cytosine, and 6-chloro-2-pyridone. Using PCM, incorrect results were obtained for isonicotinic acid and 6-chloro-2-pyridone. As shown in previous studies, local solvation structure effects, such as hydrogen bonding, are important to compute the relative conformational stabilities in protic solvents. Because RISM could capture the local solvation structure, our method successfully reproduced the experimental data.

Acknowledgement This work was supported by the Grant-in-Aid for Scientific Research (C) (Grant No. 15K05385).

Supporting Information Available Lennard-Jones parameters (Table S1); Atom labels used in Table S1 (Figure S1). This information is available free of charge via the Internet at http://pubs.acs.org.

References (1) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic: London, 1986.

13

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

(2) Hirata, F.; Sato, H.; Ten-no, S.; Kato, S. In Combined quantum mechanical and molecular mechanical methods; Gao, J., Thompson, M. A., Eds.; American Chemical Society: Washington DC, 1998. (3) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999–3094. (4) Purvis, G. D.; Bartlett, R. J. A Full Coupled-Cluster Singles and Doubles Model: The Inclusion of Disconnected triples. J. Chem. Phys. 1982, 76, 1910–1918. (5) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A Fifth-order Perturbation Comparison of Electron Correlation Theories. Chem. Phys. Lett. 1989, 157, 479–483. (6) Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular Electronic-Structure Theory; Wiley: Chichester, 2000. (7) Cammi, R. Quantum cluster theory for the polarizable continuum model. I. The CCSD level with analytical first and second derivatives. J. Chem. Phys. 2009, 131, 164104. (8) Caricato, M. CCSD-PCM: Improving upon the reference reaction field approximation at no cost. J. Chem. Phys. 2011, 135, 074113. (9) Chandler, D.; Andersen, H. C. Optimized Cluster Expansions for Classical Fluids. II. Theory of Molecular Liquids. J. Chem. Phys. 1972, 57, 1930–1937. (10) Hirata, F.; Rossky, P. J. An extended RISM equation for molecular polar fluids. Chem. Phys. Lett. 1981, 83, 329–334. (11) Sato, H.; Hirata, F. Ab Initio Study on Molecular and Thermodynamic Properties of Water: A Theoretical Prediction of pKw over a Wide Range of Temperature and Density. J. Phys. Chem. B 1999, 103, 6596–6604.

14

ACS Paragon Plus Environment

Page 14 of 19

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

Journal of Chemical Theory and Computation

(12) Sato, H.; Hirata, F. Equilibrium and Nonequilibrium solvation structure of hexaammineruthenium (II, III) in aqueous solution: Ab initio RISM-SCF study. J. Phys. Chem. A 2002, 106, 2300–2304. (13) Sato, H.; Sakaki, S. Comparison of Electronic Structure Theories for Solvated Molecules: RISM-SCF versus PCM. J. Phys. Chem. A 2004, 108, 1629–1634. (14) Yokogawa, D.; Sato, H.; Sakaki, S. New generation of the reference interaction site model self-consistent field method: Introduction of spatial electron density distribution to the solvation theory. J. Chem. Phys. 2007, 126, 244504. (15) Yokogawa, D. New fitting approach of electrostatic potential for stable quantum mechanical calculations using the reference interaction site model. Chem. Phys. Lett. 2013, 587, 113–117. (16) Yokogawa, D.; Arifin, Electrostatic Potential Charge including Spatial Electron Density Distribution (SEDD): Application to Biosystems. Bull. Chem. Soc. Jpn. 2017, 90, 831– 837. (17) Yokogawa, D. Time-dependent density functional theory (TD-DFT) coupled with reference interaction site model self-consistent field explicitly including spatial electron density distribution (RISM-SCF-SEDD). J. Chem. Phys. 2016, 145, 094101. (18) Bartlett, R. J.; Musiał, M. Coupled-cluster theory in quantum chemistry. Rev. Mod. Phys. 2007, 79, 291–352. (19) Yokogawa, D. Linear response approximation to reference interaction site model selfconsistent field explicitly including spatial electron density distribution. Free energy. J. Chem. Phys. 2013, 138, 164109. (20) Urban, M.; Noga, J.; Cole, S. J.; Bartlett, R. J. Towards a full CCSDT model for electron correlation. J. Chem. Phys. 1985, 83, 4041–4046. 15

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

(21) Yokogawa, D.; Sato, H.; Sakaki, S. Analytical energy gradient for reference interaction site model self-consistent field explicitly including spatial electron density distribution. J. Chem. Phys. 2009, 131, 214504. (22) Hayaki, S.; Kido, K.; Yokogawa, D.; Sato, H.; Sakaki, S. A theoretical analysis of a Diels-Alder reaction in ionic liquids. J. Phys. Chem. B 2009, 113, 8227–8230. (23) Iida, K.; Yokogawa, D.; Ikeda, A.; Sato, H. Carbon dioxide capture at the molecular level. Phys. Chem. Chem. Phys. 2009, 11, 8556–8559. (24) Yokogawa, D.; Ono, K.; Sato, H.; Sakaki, S. Theoretical study on aquation reaction of cis-platin complex: RISM-SCF-SEDD, a hybrid approach of accurate quantum chemical method and statistical mechanics. Dalton Trans. 2011, 40, 11125–11130. (25) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378–6396. (26) Dunning, Jr., T. H.; Peterson, K. A.; Wilson, A. K. Gaussian basis sets for use in correlated molecular calculations. X. The atoms aluminum through argon revisited. J. Chem. Phys. 2001, 114, 9244–9253. (27) Dunning, Jr., T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (28) Kendall, R. A.; Dunning, Jr., T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796–6806. (29) Papajak, E.; Leverentz, H. R.; Zheng, J.; Truhlar, D. G. Efficient Diffuse Basis Sets: cc-pVxZ+ and maug-cc-pVxZ. J. Chem. Theory Comput. 2009, 5, 1197–1202. 16

ACS Paragon Plus Environment

Page 16 of 19

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

Journal of Chemical Theory and Computation

(30) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09 Revision B.1. Gaussian Inc. Wallingford CT 2009. (31) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A. V.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; WilliamsYoung, D.; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J. J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Keith, T. A.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian 16 Revision A.03. Gaussian Inc. Wallingford CT 17

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

2016. (32) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347–1363. (33) Klopper, W.; Noga, J.; Koch, H.; Helgaker, T. Multiple Basis Sets in Calculations of Triples Corrections in Coupled-Cluster Theory. Theor. Chem. Acc. 1997, 97, 164–176. (34) Schwenke, D. W. The Extrapolation of One-electron Basis Sets in Electronic Structure Calculations: How It Should Work and How It Can be Made to Work. J. Chem. Phys. 2005, 122, 014107. (35) Koch, H.; Jørgensen, P. Coupled cluster response functions. J. Chem. Phys. 1990, 93, 3333–3344. (36) Albert, A.; Taylor, P. J. The Tautomerism of 1,2,3-Triazole in Aqueous Solution. J. Chem. Soc., Perkin Trans. 2 1989, 1903–1905. (37) Abboud, J.-L. M.; Foces-Foces, C.; Notario, R.; Trifonov, R. E.; Volovodenko, A. P.; Ostrovskii, V. A.; Alkorta, I.; Elguero, J. Basicity of N-H- and N-Methyl-1,2,3-triazoles in the Gas Phase, in Solution, and in the Solid State - An Experimental and Theoretical Study. Eur. J. Org. Chem. 2001, 2001, 3013–3024. (38) Nagy, P. I.; Alagona, G.; Ghio, C. J. Chem. Theory Comput. 2007, 3, 1249–1266. (39) Dreyfus, M.; Bensaude, O.; Dodin, G.; Dubois, J. J. Am. Chem. Soc. 1976, 98, 6338– 6349. (40) Forlani, L.; Cristoni, G.; Boga, C.; Todesco, P. E.; Vecchio, E. D.; Selva, S.; Monari, M. ARKIVOC (Gainesville, FL, U. S.) 2002, 198–215.

18

ACS Paragon Plus Environment

Page 18 of 19

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

Journal of Chemical Theory and Computation

Graphical TOC Entry CC equation

Relative stabilities

〈 μ| e-T H0solveT| Φ0〉=

0

This work PCM

K→E

C1→C3

N→Z

RISM equation 1H→2H

huv = ωu*cuv*(ωv+ρhvv)

19

0

2

4

Error (kcal/mol)

ACS Paragon Plus Environment

6