The SCC-DFTB-PIMD method to evaluate a multidimensional

Subscriber access provided by Nottingham Trent University

Reaction Mechanisms

The SCC-DFTB-PIMD method to evaluate a multidimensional quantum free energy surface for a proton transfer reaction Kento Kosugi, Hiroshi Nakano, and Hirofumi Sato J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00355 • Publication Date (Web): 16 Aug 2019 Downloaded from pubs.acs.org on August 19, 2019

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

The SCC-DFTB-PIMD method to evaluate a multidimensional quantum free energy surface for a proton transfer reaction Kento Kosugi,† Hiroshi Nakano,∗,†,‡ and Hirofumi Sato†,‡ †Department of Molecular Engineering, Kyoto University, Kyoto Daigaku Katsura, Kyoto 615-8510, Japan ‡Elements Strategy Initiative for Catalysts and Batteries, Kyoto University E-mail: [email protected]

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

Keywords: DFTB, PIMD, free energy surface, nuclear quantum effect, proton transfer Abstract The self-consistent charge density functional tight binding was combined with the path integral molecular dynamics method for the first time to evaluate the two dimensional free energy surface including nuclear quantum effects of a proton transfer reaction in a 2,4-dichlorophenol-trimethylamine complex. A statistically converged two dimensional quantum free energy surface was evaluated by using the multidimensional blue moon ensemble method. The accuracy was guaranteed by optimizing the repulsive potential between the sp3 hybridized nitrogen and hydrogen atoms in a SCC-DFTB3 parameter set for the system to reproduce high level quantum chemical calculations. The present study illustrates the usefulness of this new approach to investigate nuclear quantum effects in various realistic proton transfer reactions.

2


Page 2 of 39

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


1

Introduction

The nuclear quantum effects (NQEs) are important factors affecting the mechanisms and rates of some reactions, especially of proton transfer reactions. 1–9 A number of studies show that multidimensional energy surface is necessary in order to properly analyze the NQEs. 10–15 A representative approach for that purpose is path-integral molecular dynamics (PIMD) calculations. 16–20 By comparing the free energy surface with that from a corresponding classical calculation, the NQEs emerge as the change in the shape of the surface and decrease in the barrier height. In principle, any types of methods are applicable to get the potential energy and nuclear gradients in PIMD calculations . If the electronic Schrödinger equation is solved nonempirically, it is called as ab initio path-integral molecular dynamics (AI-PIMD) approach. 21–29 Its advantage is the applicability to any reactions and the ability to provide reliable results. The use of AI-PIMD method is however somewhat limited due to the large computational cost and the difficulty in conducting adequate configurational sampling in particular to obtain a free energy surface. An alternative approach is to combine the PIMD approach with semiempirical quantum chemical methods like AM1, PM3, and PM6, where computational efficiency is much higher than the ab initio methods. The quantum free energy profiles were calculated for some proton transfer reactions by using this approach. 30–34 The self-consistent charge density functional tight binding (SCC-DFTB) 35–37 is another powerful semiempirical method. The method is computationally efficient similar to other semiempirical methods. It also gives quantitative results comparable to or better than not only other semiempirical methods but also DFT calculations with small basis sets by using its transferable parameter sets. 38 Although there are some difficult systems to describe quantitatively using a general-purpose parameter set, 37–43 SCC-DFTB method enables us to sample enough configurations with much lower computational cost to obtain statistically converged multidimensional free energy surface. The accuracy can be further improved by optimizing a part of the parameter set, especially the pair repulsive potentials, for the specific 3



system being studied. In this study, we will present SCC-DFTB-PIMD approach (PIMD simulations using SCCDFTB to obtain energy and nuclear gradients) to calculate a two dimensional quantum free energy surface for a proton transfer reaction. The free energy surface was evaluated by using the multidimensional blue moon ensemble method. 44–46 As far as we know, this is the first attempt to combine SCC-DFTB with PIMD, and also with the multidimensional blue-moon ensemble method to calculate a quantum free energy surface for a proton transfer reaction. We also found that SCC-DFTB3, the newest version of SCC-DFTB at present, provides quantitative results in good agreement with those of high level quantum chemical methods by optimizing a pair repulsive potential. Here we examined a proton transfer reaction in a 2,4-dichlorophenol-trimethylamine complex in the gas phase. Although there have been many important studies of NQEs for phenol-amine complexes 47,48 by using empirical potential functions, 49–52 only a limited number of quantum chemical calculations are reported to quantitatively study a realistic molecular system. 53,54 The reason is probably the system is somewhat more complex compared to another well known system such as an intramolecular proton transfer in malonaldehyde, where the size is smaller and the barrier height is significantly lower than that of the present system. In order to obtain a statistically converged two-dimensional free energy surface, the multidimensional blue moon ensemble method was employed. The NQEs were analyzed by focusing on the extent of expansion of the transferring proton ring polymer.

4


Page 4 of 39

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

Methods

2.1

SCC-DFTB-PIMD multidimensional blue-moon ensemble approach

The equilibrium statistical quantities can be evaluated quantum mechanically from the partition function

( )P ˆ ˆ = lim ZP (β), Z(β) = Tre−β H = Tr e−β H/P P →∞

(1)

ˆ = H(ˆ ˆ rN , p ˆ N ) is the quantum mechanical Hamiltonian of the system where β = 1/kB T and H composed of N quantum nuclei. The discretized partition function ZP (β) is given by ∫ ZP (β) =

drN P dpN P e−βHP (r

N P ,pN P )

(2)

with the classical mechanical Hamiltonian HP (rN P , pN P ) of the system composed of N ring polymers of P beads expressed as [ N ] P N ∑ p2i,k ∑ ∑ 1 1 HP (rN P , pN P ) = + mi ωP2 (ri,k+1 − ri,k )2 + EDFTB ({ri,k }N i=1 ) . ′ 2m 2 P i i=1 k=1 i=1

(3)

Here ri,P +1 = ri,1 , pi,P +1 = pi,1 , and {ri,k }N i=1 ≡ r1,k , r2,k , · · · , rN,k . mi is the ith particle √ mass, m′i = mi P/(2πℏ)2 , and ωP = P /βℏ is the frequency of the harmonic spring interconnecting adjacent polymer beads. 16–20 EDFTB ({ri,k }N i=1 ) is the potential energy evaluated with SCC-DFTB3 for the molecular system composed of N nuclei, whose coordinates are {ri,k }N i=1 . The phase space points for N ring polymers of P beads are sampled from classical MD calculations using the Hamiltonian in Eq. (3). The potential energy EDFTB ({ri,k }N i=1 ) is thus needed to be evaluated P times at each MD time step. The SCC-DFTB3 potential energy of a molecule consisting of N nuclei can be written

5



Page 6 of 39

as 37

EDFTB (rN ) = E 0 (rN ) + E 2nd (rN ) + E 3rd (rN ) + E rep (rN ) occ. ∑ ∑ 1∑ ∗ 0 = ni Cµi Cνi Hµν + γIJ (|rI − rJ |)∆qI ∆qJ 2 µ,ν i I,J 1∑ 1 ∑ rep + ΓIJ (|rI − rJ |)∆qI2 ∆qJ + V (|rI − rJ |). 3 I,J 2 I,J IJ 0 Here Hµν =

⟨

(4)

⟩ ˆ 0 ϕµ H ϕν , Sµν = ⟨ϕµ |ϕν ⟩ with {ϕµ } being valence pseudoatomic orbitals,

rep γIJ , ΓIJ , and VIJ are functions of distance between atom pairs precalculated and tabulated, ∑ ∑ and these are considered as parameters. ∆qI = µ∈I ν Pµν Sµν − ZI is the Mulliken charge ∑ ∗ on Ith nuclei (taking care of a sign convention), where Pµν = occ i ni Cµi Cνi is the density

matrix, Cµi is ith molecular orbital coefficient on µth pseudoatomic orbital, and ZI is the core charge on Ith nuclei. Minimizing Eq. (4) with respect to {Cµi } gives a secular equation to determine {Cµi } in a self-consistent manner. The reference energy variations for a dimer ref EIJ (|rI − rJ |) as a function of the distance |rI − rJ | are reproduced by determining the

repulsive potential as follows { 0 } rep ref 2nd 3rd VIJ (|rI − rJ |) = EIJ (|rI − rJ |) − EIJ (|rI − rJ |) + EIJ (|rI − rJ |) + EIJ (|rI − rJ |) . (5)

We focus on modifying only this respulsive potential between nitrogen and hydrogen atoms in an existing parameter set without changing the other parameters. The two dimensional free energy surface along reaction coordinates ξ1 and ξ2 are evaluated by integrating the mean forces as follows ∫ F (ξ1 , ξ2 ) −

F (ξ10 , ξ20 )

ξ1

= ξ10

6

∂F ′ dξ + ∂ξ1′ 1

∫


ξ2

ξ20

∂F ′ dξ . ∂ξ2′ 2

(6)

Page 7 of 39 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 mean forces are calculated from 44–46 ⟨ −1/2 ⟩ |Z| (kB T Gα /2 − λα ) ξ′ ∂F = ∂ξα′ ⟨|Z|−1/2 ⟩ξ′

(α = 1, 2)

(7)

where ⟨X⟩ξ′ indicates that the quantity X is evaluated in a simulation constraining the reaction coordinates as (ξ1 , ξ2 ) = (ξ1′ , ξ2′ ). λα is the Lagrange multiplier to constrain αth reaction coordinate represented as |rNH (t + ∆t)| − |rOH (t + ∆t)| − ξ1′ , ∆ˆrO,N (t) · ˆrNH (t + ∆t) + ∆ˆrN,O (t) · ˆrOH (t + ∆t) ξ2′2 − |rON (t + ∆t)|2 = −1 rNO (t) 2(m−1 O + mN )rON (t + ∆t) · ˆ

λ1 =

(8)

λ2

(9)

obtained from the RATTLE algorithm, where rIJ = rI − rJ , ˆrIJ = rIJ /|rIJ |, and ∆ˆrI,J = −1 m−1 rIH − (m−1 rJH with {I, J} = {O, N, H}. The 2 × 2 matrix Z and the vector G H ˆ J + mH )ˆ

are respectively defined as follows 46

Zαβ

N ∑ 1 ∂ξα ∂ξβ = · mi ∂ri ∂ri i=1

Gα =

N ∑ i,j

(α, β = 1, 2)

2 ∑ 2 ∂ξδ ∂ 2 ξϵ ∂ξγ −1 (Z −1 )δϵ · · (Z )γα mi mj γ,δ,ϵ=1 ∂ri ∂ri ∂rj ∂rj

(10)

(α = 1, 2)

(11)

In the present case, G1 = G2 = 0 since ∂ 2 ξα /∂ri ∂rj = 0 (α = 1, 2), and 1 1 2 + + (1 − ˆrNH · ˆrOH ) mO mN mH 1 1 ˆrOH · ˆrON + ˆrNH · ˆrNO = Z21 = − mO mN 1 1 = + . mO mN

Z11 =

(12)

Z12

(13)

Z22

(14)

Note that the mass of the centroid in the PIMD calculations is equal to the actual atom mass as explained below. The Cartesian coordinates of the ring-polymer beads are transformed to the normal-mode 7



Page 8 of 39

coordinates and evolved in time to accelerate the phase space sampling. The normal modes can be obtained by using the transformation matrix C as follows 19,20

ui,k

P 1 ∑ = Ck,l ri,l P l=1

(15)

ri,k

P 1 ∑ T = (C )k,l ui,l . P l=1

(16)

The Hamiltonian in Eq. (3) is then rewritten in the normal-mode representation with the corresponding momenta as [ N ] P N ∑ ∑ p2i,k ∑ 1 1 HP (uN P , pN P ) = + mi,k ωP2 u2i,k + EDFTB ({ri,k (uN )}N i=1 ) . ′ 2m 2 P i,k i=1 k=1 i=1

(17)

The transformation matrix is defined for even P as

Cj,k =

 √    1/P , k=1      √    2/P cos(2πjk/P ), 2 ≤ k ≤ P/2 √    1/P (−1)j ,      √    2/P sin(2πjk/P ),

(18)

k = P/2 + 1 P/2 + 2 ≤ k ≤ P

and for odd P as

Cj,k

 √    1/P , k=1     √ = 2/P cos(2πjk/P ), 2 ≤ k ≤ (P + 1)/2     √    2/P sin(2πjk/P ), (P + 1)/2 + 1 ≤ k ≤ P.

(19)

The masses in Eq. (17) are given with the eigenvalues of the transformation matrix C, λ2k−1 = λ2k = 2[1 − cos{2π(k − 1)/P }] like mi,k = mi λk , m′i,1 = mi , and m′i,k = mi,k (k = ∑ 2, 3, · · · , P ). Note that the mass mi,1 = mi associated with the centroid ui,1 = Pk=1 ri,k /P

8


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


is equal to that of the actual mass of ith nucleus, which is used in the blue-moon ensemble calculations in Eqs.(12), (13), and (14).

3

Computational Details

The optimized reactant structure of the 2,4-dichlorophenol-trimethylamine dimer handled in this study is displayed in Fig. 1(a). The following reaction coordinates were employed

The centroid ui,1 =

∑P k=1

ξ1 = r = |rN − rH | − |rO − rH |,

(20)

ξ2 = R = |rN − rO |.

(21)

ri,k /P was used for these reaction coordinates in the SCC-DFTB-

PIMD blue-moon ensemble calculations. The 3ob parameter set 38 was used in the SCC-DFTB3 potential energy and nuclear rep gradient calculations. The N-H pair repulsive potential VNH was reoptimized so that it

reproduces the potential energy profile between a trimethylamine N(CH3 )3 and a proton calculated at the B3LYP/6-31G(d,p) level 55–57 based on Eq. (5) from r(N-H) = 0.95 ˚ A to r(N-H) = 2.33 ˚ A. SCC-DFTB-MD (classical MD simluations using SCC-DFTB3 to obtain energy and nuclear gradients) and SCC-DFTB-PIMD calculations were performed in the gas phase by using the velocity Verlet and RATTLE algorithms 58 with 0.5 fs timestep in the NVT ensemble at 249.0 K as in the previous works. 49–54 The temperature was controlled by using the massive Nosé-Hoover chain thermostats. 59 The MD trajectories were propagated using the side scheme, 60 in which the phase space propagator eL∆t is evaluated as eL∆t ∼ eLT ∆t/2 eLp ∆t/2 eLr ∆t eLp ∆t/2 eLT ∆t/2 ,

(22)

where LT , Lp , and Lr are Liouville operators for the thermostat, momenta, and positions, 9



respectively. It should be noted that a previous work showed the importance of using a smaller timestep to obtain converged kinetic and potential energies in PIMD simulations. 60 We thus also used 0.1 fs timestep to calculate the one dimensional free energy profiles. At the same time, the convergence of the barrier height with respect to timestep was examined along the profiles with various R values in reference to another previous work. 61 The results convince us that the convergence is safely achieved with 0.5 fs timestep in this system (see Appendix). In the PIMD calculations, the Cartesian coordinates of the beads were transformed to the normal mode coordinates as shown in the previous section and evolved in time to accelerate the phase space sampling. 19,20 The reaction coordinates were discretized into grid points of 0.10 ˚ A width unless otherwise noted. At each grid point, 10 ps production run was conducted to obtain the mean force for the free enegy calculation by thermodynamic integration. The statistical convergence was carefully checked by comparing the minimum free energy curves calculated with longer production runs up to 50 ps as shown in Appendix. The number of beads (the imaginary time slices) P was set to 16. The convergence with respect to the beads number was examined as illustrated in Appendix, and the error in the NQEs on the free energy surface was estimated to be within 0.2 kcal/mol. All of the calculations were performed with a modified version of the GAMESS program package. 62

4 4.1

Results and Discussion Performance of the optimized N-H pair repulsive potential

We first examined the performance of the SCC-DFTB3 calculation of the potential energy profile for the present proton transfer reaction. Fig. 2 illustrates the profiles along the reaction coordinate r in Eq. (20) calculated by using the 3ob parameter set with different N-H pair repulsive potentials; the original one, 38 a modified one taken from previous works (NHmod), 37–40 and that optimized for the present system (NHopt). The results from MP2 and B3LYP calculations with the 6-31G(d,p) basis set are also shown for comparison, which 10


Page 10 of 39

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


give almost the same profiles and overlap each other. All the degrees of freedom other than the reaction coordinate were optimized in the calculations. The calculation with the original N-H pair repulsive potential gives the energy difference between the product side (r ∼ −1.0 ˚ A) and the reactant side (r ∼ 1.0 ˚ A) about 9 kcal/mol higher than those at the reference MP2 and B3LYP levels. This result is consistent with previous reports where considerable errors can be observed when a molecule containing sp3 hybridized nitrogen is treated. 39–43 The energy difference is considerably improved by using the NHmod repulsive potential. However, nonnegligible error remains in the reactant region (0.0 ˚ A ≲ r ≲ 1.0 ˚ A), which can have a large impact on the analysis of the reaction. On the other hand, the new NHopt repulsive potential optimized for the present system gives the profile in almost perfect agreement with the reference profiles. The deviation is at most 0.2 kcal/mol. The N-H repulsive potential shown in Fig. 3. NHmod makes a N-H bond more stable than the original one as seen as the deeper minimum around r ∼ 1.2 ˚ A. NHopt has a broader minimum compared to NHmod and gives a more attractive interaction between N and H atoms in somewhat longer range, resulting in the improvement in the potential energy profile in the reactant region for this specific system. We also compared the two dimensional potential energy surface for 2.4 ˚ A ≤ R ≤ 3.1 ˚ A and −1.0 ˚ A ≤ r ≤ 1.0 ˚ A calculated at the SCC-DFTB3 with NHopt and the B3LYP/6-31G(d,p) levels (the figures are not shown since they look similar to the free energy surfaces in Fig. 4). The largest difference is found to be 1.5 kcal/mol in the region of r ∼ 0.0 ˚ A and R ≥ 2.8 ˚ A, and they are overall in quantitative agreement with each other. We also validated NHopt by comparing the optimized reactant structures. Table 1 lists two characteristic distances (r = |rN −rH |−|rO −rH | and R = |rN −rO |) and the angle formed by the rHO and rHN vectors in the optimized geometry. Although the angles are similar to each other, the distances are considerably different among the SCC-DFTB3 calculations with three different N-H repulsive potentials. The original one overestimates r and R by about 0.3 ˚ A compared to the B3LYP results. The use of NHmod produces almost the same

11



distances as the original one. NHopt significantly improves the optimized structure, and the deviations in the distances compared to the reference B3LYP results are less than 0.03 ˚ A. These results indicate the N-H pair repulsive potential optimized for the present system accurately gives not only the potential energies but also the structures, and can be used for quantitative analysis on the NQEs in the proton transfer reaction.

4.2

Classical and quantum free energy surfaces

The two dimensional free energy surfaces were calculated with the SCC-DFTB-MD and SCC-DFTB-PIMD methods as shown in Fig. 4 (a) and (b), respectively. The free energy value at the minimum was set to zero in each surface plot. The bottom around the lower right corner in (a) and (b) corresponds to the stable reactant. The free energy increases as the reaction coordinate r decreases toward the product at the upper left corner. Although the overall shape of the surfaces look similar, nonnegligible NQEs are found. The difference between (a) and (b) is also shown in (c), which corresponds to the NQEs on the free energy surface. Because of the large quantum effects when R is large and r ∼ 0.0 ˚ A, the free energy barrier disappears along the path with R fixed at 2.8 ˚ A as shown in Fig. 5. Fig. 4 (a) and (b) also depict with broken lines the minimum free energy paths obtained from the onedimensional free energy calculations along r. The NQEs modify slightly the minimum free energy path and considerably the position of the free energy minimum as rmin = 1.0 ˚ A in the SCC-DFTB-MD calculation while rmin = 0.88 ˚ A in the SCC-DFTB-PIMD calculation. The NQEs also affect the free energy values along the path as shown in Fig. 6, in which the difference between them is also depicted. Although the classical and quantum free energy calculations give almost the same free energy difference between the product side and the reactant side, we can find difference larger than 2 kcal/mol in the intermediate region (r ∼ 0.0 ˚ A), where the distances r(N-H) and r(O-H) are about 1.25 ˚ A in both the SCC-DFTB-MD and the SCC-DFTB-PIMD calculations. The degree of NQEs on the free energy profile along the reaction path is comparable 12


Page 12 of 39

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


to that reported in previous calculations on a corresponding model system in solution. 51,52 Note that the potential (or free) energy in the product side is much smaller in the previous calculations and there exists a double well in the profile. Namely, there is a large difference in the shape of the potential energy surface based on which the quantum free energy calculations were performed. Our calculations used a quantum chemical method to evaluate the potential energy and free energy of a real system, though the solvation effects are not included in the present work.

4.3

Expansion and contraction of the proton ring polymer

It is interesting to consider the reason why the free energy change due to the NQEs becomes larger on a path with larger R = |rN − rO | as seen in Fig. 4 (c). Since the NQEs on the free energy profile have relation to the variation in the degree of expansion of the proton ring polymer along the path, we evaluated the average distance between the proton beads and their centroid defined as (see also Fig. 7(a)) P P 1 ∑ 1 ∑ C |r | = |rH,k − uH,1 |, ∆|rH | = P k=1 k P k=1

(23)

where uH,1 is the centroid of the proton ring polymer defined under Eq. (21). This quantity becomes large when the ring polymer is expanded, meaning that ∆|rH | corresponds to a quantity to represent the delocalization of a quantum particle. 63–65 Fig. 8 plots ∆|rH | along the paths with R = 2.3, 2.7, and 3.1 ˚ A. The extent of the expansion of the proton ring polymer remains small and almost the same over the path with R = 2.3 ˚ A. On the other hand, ∆|rH | considerably changes and the proton ring polymer becomes more expanded with R = 2.7 ˚ A and 3.1 ˚ A as |r| departs from 1.0 ˚ A, where the O-H bond or the N-H bond is broken. These results indicate the NQEs are rather small along the path with R = 2.3 ˚ A, but they become larger as R increases like R = 2.7 ˚ A and 3.1 ˚ A. This can be verified from Fig. 4 (c): The free energy profile along the path with R = 2.3 ˚ A

13



exhibits small NQEs by about 1 kcal/mol. Larger NQEs by more than 3 kcal/mol emerge along the path with R = 2.7 ˚ A and by more than 5 kcal/mol along the path with R = 3.1 ˚ A. Since R along the minimum free energy path varies like 2.5 ˚ A ≲ R ≲ 2.9 ˚ A as seen in Fig. 4 (a), (b), the variation of ∆|rH | and the NQEs on the path are expected to be similar to those along the path with R = 2.7 ˚ A fixed, as verified from Fig. 6 for the NQEs and Fig. 9 for the extent of ring-polymer expansion with broken line. The expansion of the ring polymer should be related to the force acting on the beads. We thus evaluated the force ∂EDFTB ({ri,k }N i=1 )/∂rH,k projected onto the vector from each bead to the centroid as illustrated in Fig. 7(b). The ring polymer becomes contracted as the projected force becomes large. Fig. 9 illustrates the projected force averaged over all beads of the proton ring polymer and also over the trajectory at each point on the minimum free energy path. A clear correlation can be found between the projected force and the extent of the expansion (contraction) of the ring polymer. The force is strong at the reactant side and the product side of |r| ∼ 1.0 ˚ A. This is due to the bonding force between the proton and the oxygen or the nitrogen atom to move the beads to a stable equilibrium position where the centroid tends to exist at the reactant and the product regions. The bonding force becomes weak at r ∼ 0.0 ˚ A, where the O-H bond or the N-H bond is at least partially broken. This trend emerges more clearly on the path with R = 3.1 ˚ A, where both of the O-H and N-H bonds are almost completely broken at around r = 0.0 ˚ A. On the other hand, the bonds remain (partially) over the path with R = 2.3 ˚ A and then the variation in the contraction force is small. These analysis suggest that the degree of the NQEs on the free energy profile can be estimated from the variation in the strength of the force acting on the ring-polymer beads to contract to the centroid position. The contraction force mainly comes from the O-H and N-H bonds in the case of this reaction in the gas phase. The variation of the contraction force should also be related to the variation in the zeropoint energies. Note that the zero-point energies in the directions orthogonal to the reaction coordinate contribute to all over the reaction profile, even when a minimum is absent at the

14


Page 14 of 39

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


product region and the profile is not in accordance with a conventional double well picture. The zero-point energies are large at the reactant and product regions because the bonds produce strong contraction force not only in the bonding direction but also in some of the orthogonal directions. The zero-point energies become smaller as r approaches zero as the bonding force becomes smaller. In order to estimate the zero-point energy contributions to the NQEs along the minimum free energy path, the following projected Hessian was used

K P (ξ1 ) = (1 − P (ξ1 ))K(ξ1 )(1 − P (ξ1 ))

(24)

to get the vibrational modes orthogonal to the reaction coordinate and their zero-point energies. Here K(ξ1 ) is the Hessian and P (ξ1 ) is the matrix to project out the components parallel to the reaction coordinate. 66 The resulting zero-point energies were 128.7, 124.0, and 128.0 kcal/mol at r = −1.0, 0.0, and 1.0 ˚ A, respectively. These reduce the (classical) free energy difference between r = 0.0 ˚ A and r = 1.0 ˚ A by 4.0 kcal/mol, and between r = −1.0 ˚ A and r = 1.0 ˚ A by -0.7 kcal/mol, which slightly overestimates the NQEs but comparable to those (2.5 and -0.1 kcal/mol, respectively) shown in Fig. 6 (b).

5

Concluding remarks

In the present study, we combined SCC-DFTB3 with PIMD for the first time to evaluate the two dimensional free energy surface for the proton transfer reaction in a 2,4-dichlorophenoltrimethylamine complex in gas phase. The statistically converged free energy surface was calculated by using the multidimensional blue-moon ensemble method. The accuracy of the free energies was ensured by optimizing the N-H pair repulsive potential in the SCC-DFTB3 energy and nuclear gradient calculations for this specific system. The use of this optimized repulsive potential reproduced the reference results evaluated at the B3LYP and MP2 levels quite well. The NQEs on the reaction paths are small when the distance R between the oxygen and nitrogen atoms is small. The reaction path with large R exhibits large NQEs. 15



The NQEs are also observed on the minimum free energy path, though it does not have a double well as in the cases discussed in most of the studies on the NQEs on reactions. In order to analyze these results, we evaluated the degree of expansion of the transferring proton ring polymer and its variation along the reaction path. It was found that the extent of the ring polymer expansion considerably changes along the reaction path with large R, and the opposite trend was observed along the path with small r. A clear correlation was also observed between this degree of the ring polymer expansion and the strength of the contraction force that drives the ring polymer beads to the centroid, which mainly originates from the O-H or N-H chemical bonds in the present case. These factors, especially about the contraction force, were interpreted with relevance to the zero-point energies. These analysis provides a way to understand how the NQEs appear on a reaction path. Because of the good accuracy and the low computational cost, the present SCC-DFTBPIMD approach can serve as a powerful tool to investigate reactions where the nuclear quantum effects play an important role in detail, so long as some repulsive potentials used in the SCC-DFTB3 calculations are appropriately optimized for the specific system if necessary. An extension of the present approach to a SCC-DFTB3/MM-PIMD to deal with condensed-phase reactions is straightforward and currently in progress in our laboratory, and the corresponding calculations in solution will be presented in the near future.

Acknowledgement The authors thank Dr. Yoshio Nishimoto for helpful discussion on the SCC-DFTB3 method. This work was financially supported by Grant-in-Aid for Young Scientists (B) (18K14179) and by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) Japan. A part of this work was performed under a management of Elements Strategy Initiative for Catalysts and Batteries (ESICB).

16


Page 16 of 39

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


6

Appendix: Convergence of the calculated free energies

The convergence with respect to the number of beads in the SCC-DFTB-PIMD free energies was examined by varying the number of beads from P = 1 to 64 as illustrated in Fig. 10. The calculation with P = 16 used in the present study gives the free energy profile in good agreement with that with P = 64 showing at most 0.2 kcal/mol deviation over the profile. We also carefully checked the statistical convergence of the free energy values. Fig. 11 displays the SCC-DFTB-PIMD (P = 16) free energy profiles obtained from 10, 20, 30, 40, and 50 ps trajectories. The standard deviations are within 0.05 kcal/mol throughout the profile. The convergence of the free energy profile and the barrier height were examined with respect to the timestep used in the SCC-DFTB-PIMD calculations (P = 16, 10 ps trajectories). Fig. 12 shows the free energy profiles calculated with various timesteps at R = 2.8 ˚ A, where the free energy barrier appears around r = −0.2 ˚ A in the SCC-DFTB-MD calculation as shown in Fig. 5. The corresponding barrier heights were extracted into Fig. 13, which are defined as the difference between the free energies at r = −0.2 ˚ A and r = 0.9 ˚ A. The results at various R values indicate that the difference is within 0.2 kcal/mol between the barrier heights calculated with 0.5 fs timestep and with 0.1 fs timestep.

17



References (1) Al Taylor, C.; El-Bayoumi, M. A.; Kasha, M. Excited-state two-proton tautomerism in hydrogen-bonded N-heterocyclic base pairs. Proc. Natl. Acad. Sci. USA 1969, 63, 253–260. (2) Okamura, M.; Feher, G. Proton transfer in reaction centers from photosynthetic bacteria. Ann. Rev. Biochem. 1992, 61, 861–896. (3) Marx, D. Proton transfer 200 years after von Grotthuss: Insights from ab initio simulations. Chem. Phys. Chem. 2006, 7, 1848–1870. (4) Umena, Y.; Kawakami, K.; Shen, J.-R.; Kamiya, N. Crystal structure of oxygenevolving photosystem II at a resolution of 1.9 ˚ A. Nature 2011, 473, 55. (5) Weinberg, D. R.; Gagliardi, C. J.; Hull, J. F.; Murphy, C. F.; Kent, C. A.; Westlake, B. C.; Paul, A.; Ess, D. H.; McCafferty, D. G.; Meyer, T. J. Proton-coupled electron transfer. Chem. Rev. 2012, 112, 4016–4093. (6) Saito, K.; Rutherford, A. W.; Ishikita, H. Mechanism of proton-coupled quinone reduction in Photosystem II. Proc. Nat. Acad. Sci. USA 2013, 110, 954–959. (7) Inagaki, T.; Yamamoto, T. Critical role of deep hydrogen tunneling to accelerate the antioxidant reaction of ubiquinol and vitamin E. J. Phys. Chem. B 2014, 118, 937–950. (8) Meisner, J.; Kästner, J. Atom tunneling in chemistry. Angew. Chem. Int. Ed. 2016, 55, 5400–5413. (9) Sakaushi, K.; Lyalin, A.; Taketsugu, T.; Uosaki, K. Quantum-to-Classical Transition of Proton Transfer in Potential-Induced Dioxygen Reduction. Phys. Rev. Lett. 2018, 121, 236001. (10) Babamov, V.; Marcus, R.-A. Dynamics of hydrogen atom and proton transfer reactions. Symmetric case. J. Chem. Phys. 1981, 74, 1790–1798. 18


Page 18 of 39

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


(11) Truhlar, D. G.; Garrett, B. C. Dynamical bottlenecks and semiclassical tunneling paths for chemical reactions. J. Chem. Phys. 1987, 84, 365–369. (12) Shida, N.; Barbara, P. F.; Almlöf, J. E. A theoretical study of multidimensional nuclear tunneling in malonaldehyde. J. Chem. Phys. 1989, 91, 4061–4072. (13) Shida, N.; Almlof, J.; Barbara, P. F. Tunneling paths in intramolecular proton transfer. J. Phys. Chem. 1991, 95, 10457–10464. (14) Tautermann, C. S.; Voegele, A. F.; Loerting, T.; Liedl, K. R. The optimal tunneling path for the proton transfer in malonaldehyde. J. Chem. Phys. 2002, 117, 1962–1966. (15) Meuwly, M.; M¨ uller, A.; Leutwyler, S. Energetics, dynamics and infrared spectra of the DNA base-pair analogue 2-pyridone· 2-hydroxypyridine. Phys. Chem. Chem. Phys. 2003, 5, 2663–2672. (16) Chandler, D.; Wolynes, P. G. Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids. J. Chem. Phys. 1981, 74, 4078–4095. (17) Parrinello, M.; Rahman, A. Study of an F center in molten KCl. J. Chem. Phys. 1984, 80, 860–867. (18) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. Efficient molecular dynamics and hybrid Monte Carlo algorithms for path integrals. J. Chem. Phys. 1993, 99, 2796–2808. (19) Cao, J.; Martyna, G. Adiabatic path integral molecular dynamics methods. II. Algorithms. J. Chem. Phys. 1996, 104, 2028–2035. (20) Tuckerman, M. E. Statistical mechanics: theory and molecular simulation; Oxford university press, 2010. (21) Marx, D.; Parrinello, M. Ab initio path integral molecular dynamics: Basic ideas. J. Chem. Phys. 1996, 104, 4077–4082. 19



(22) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M. Efficient and general algorithms for path integral Car–Parrinello molecular dynamics. J. Chem. Phys. 1996, 104, 5579–5588. (23) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M. On the quantum nature of the shared proton in hydrogen bonds. Science 1997, 275, 817–820. (24) Miura, S.; Tuckerman, M. E.; Klein, M. L. An ab initio path integral molecular dynamics study of double proton transfer in the formic acid dimer. J. Chem. Phys. 1998, 109, 5290–5299. (25) Tuckerman, M. E.; Marx, D. Heavy-atom skeleton quantization and proton tunneling in intermediate-barrier hydrogen bonds. Phys. Rev. Lett. 2001, 86, 4946. (26) Shiga, M.; Tachikawa, M.; Miura, S. A unified scheme for ab initio molecular orbital theory and path integral molecular dynamics. J. Chem. Phys. 2001, 115, 9149–9159. (27) Kaczmarek, A.; Shiga, M.; Marx, D. Quantum effects on vibrational and electronic spectra of hydrazine studied by on-the-fly ab initio ring polymer molecular dynamics. J. Phys. Chem. A 2009, 113, 1985–1994. (28) Ogata, Y.; Daido, M.; Kawashima, Y.; Tachikawa, M. Nuclear quantum effects on protonated lysine with an asymmetric low barrier hydrogen bond: an ab initio path integral molecular dynamics study. RSC Advances 2013, 3, 25252–25257. (29) Machida, M.; Kato, K.; Shiga, M. Nuclear quantum effects of light and heavy water studied by all-electron first principles path integral simulations. J. Chem. Phys. 2018, 148, 102324. (30) Major, D. T.; Garcia-Viloca, M.; Gao, J. Path integral simulations of proton transfer reactions in aqueous solution using combined QM/MM potentials. J. Chem. Theory Comput. 2006, 2, 236–245. 20


Page 20 of 39

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


(31) Zimmermann, T.; Van´ıˇcek, J. Path integral evaluation of equilibrium isotope effects. J. Chem. Phys. 2009, 131, 024111. (32) Yoshikawa, T.; Sugawara, S.; Takayanagi, T.; Shiga, M.; Tachikawa, M. Theoretical study on the mechanism of double proton transfer in porphycene by path-integral molecular dynamics simulations. Chem. Phys. Lett. 2010, 496, 14–19. (33) Suzuki, K.; Kayanuma, M.; Tachikawa, M.; Ogawa, H.; Nishihara, H.; Kyotani, T.; Nagashima, U. Path integral molecular dynamics for hydrogen adsorption site of zeolitetemplated carbon with semi-empirical PM3 potential. Comput. Theor. Chem. 2011, 975, 128–133. (34) Sugawara, S.; Yoshikawa, T.; Takayanagi, T.; Shiga, M.; Tachikawa, M. Quantum proton transfer in hydrated sulfuric acid clusters: A perspective from semiempirical path integral simulations. J. Phys. Chem. A 2011, 115, 11486–11494. (35) Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R. Construction of tight-binding-like potentials on the basis of density-functional theory: Application to carbon. Phys. Rev. B 1995, 51, 12947. (36) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys. Rev. B 1998, 58, 7260. (37) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: extension of the self-consistent-charge densityfunctional tight-binding method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7, 931–948. (38) Gaus, M.; Goez, A.; Elstner, M. Parametrization and benchmark of DFTB3 for organic molecules. J. Chem. Theory Comput. 2013, 9, 338–354.

21



(39) Bondar, A.-N.; Fischer, S.; Smith, J. C.; Elstner, M.; Suhai, S. Key role of electrostatic interactions in bacteriorhodopsin proton transfer. J. Am. Chem. Soc. 2004, 126, 14668– 14677. (40) Elstner, M. The SCC-DFTB method and its application to biological systems. Theoret. Chem. Acc. 2006, 116, 316–325. (41) Yang, Y.; Yu, H.; York, D.; Cui, Q.; Elstner, M. Extension of the self-consistent-charge density-functional tight-binding method: third-order expansion of the density functional theory total energy and introduction of a modified effective coulomb interaction. J. Phys. Chem. A 2007, 111, 10861–10873. (42) Vuong, V. Q.; Akkarapattiakal Kuriappan, J.; Kubillus, M.; Kranz, J. J.; Mast, T.; Niehaus, T. A.; Irle, S.; Elstner, M. Parametrization and Benchmark of Long-Range Corrected DFTB2 for Organic Molecules. J. Chem. Theory Comput. 2017, 14, 115–125. (43) Gruden, M.; Andjeklović, L.; Jissy, A. K.; Stepanović, S.; Zlatar, M.; Cui, Q.; Elstner, M. Benchmarking density functional tight binding models for barrier heights and reaction energetics of organic molecules. J. Comput. Chem. 2017, 38, 2171–2185. (44) Carter, E.; Ciccotti, G.; Hynes, J. T.; Kapral, R. Constrained reaction coordinate dynamics for the simulation of rare events. Chem. Phys. Lett. 1989, 156, 472–477. (45) Sprik, M.; Ciccotti, G. Free energy from constrained molecular dynamics. J. Chem. Phys. 1998, 109, 7737–7744. (46) Sergi, A.; Ciccotti, G.; Falconi, M.; Desideri, A.; Ferrario, M. Effective binding force calculation in a dimeric protein by molecular dynamics simulation. J. Chem. Phys. 2002, 116, 6329–6338. (47) Ilczyszyn, M.; Ratajczak, H.; Ladd, J. A. Reversible proton transfer phenomenon

22


Page 22 of 39

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


in the 2, 4-dichlorophenol-triethylamine hydrogen-bonded complex studied by lowtemperature 1H NMR spectroscopy. Chem. Phys. Lett. 1988, 153, 385–388. (48) Ilczyszyn, M.; Ratajczak, H.; Ladd, J. A. Reversible proton transfer in the 2, 3, 5, 6-tetrachlorophenol-N, N-dimethylaniline hydrogen-bonded complex studied by lowtemperature 1H NMR spectroscopy. J. Mol. Struct. 1989, 198, 499–503. (49) Azzouz, H.; Borgis, D. A quantum molecular-dynamics study of proton-transfer reactions along asymmetrical H bonds in solution. J. Chem. Phys. 1993, 98, 7361–7374. (50) Hammes-Schiffer, S.; Tully, J. C. Proton transfer in solution: Molecular dynamics with quantum transitions. J. Chem. Phys. 1994, 101, 4657–4667. (51) Yamamoto, T.; Miller, W. H. Path integral evaluation of the quantum instanton rate constant for proton transfer in a polar solvent. J. Chem. Phys. 2005, 122, 044106. (52) Collepardo-Guevara, R.; Craig, I. R.; Manolopoulos, D. E. Proton transfer in a polar solvent from ring polymer reaction rate theory. J. Chem. Phys. 2008, 128, 144502. (53) Aono, S.; Kato, S. Proton transfer in phenol–amine complexes: Phenol electronic effects on free energy profile in solution. J. Comput. Chem. 2010, 31, 2924–2931. (54) Aono, S.; Yamamoto, T.; Kato, S. Solution reaction space Hamiltonian based on an electrostatic potential representation of solvent dynamics. J. Chem. Phys. 2011, 134, 144108. (55) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098. (56) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785. (57) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. 23



(58) Anderson, H. C. Rattle: A velocity version of the shake algorithm for molecular dynamics calculations. J. Comp. Phys 1983, 52, 24–34. (59) Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé–Hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 1992, 97, 2635–2643. (60) Zhang, Z.; Liu, X.; Chen, Z.; Zheng, H.; Yan, K.; Liu, J. A unified thermostat scheme for efficient configurational sampling for classical/quantum canonical ensembles via molecular dynamics. J. Chem. Phys. 2017, 147, 034109. (61) Zhang, Z.; Liu, X.; Yan, K.; Tuckerman, M. E.; Liu, J. A Unified Efficient Thermostat Schemefor the Canonical Ensemblewith Holonomic or Isokinetic Constraints via Molecular Dynamics. J. Phys. Chem. A 2019, (62) 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. (63) Landman, U.; Scharf, D.; Jortner, J. Electron localization in alkali-halide clusters. Phys. Rev. Lett. 1985, 54, 1860. (64) Barnett, R.; Landman, U.; Cleveland, C.; Jortner, J. Electron localization in water clusters. II. Surface and internal states. J. Chem. Phys. 1988, 88, 4429–4447. (65) Benoit, M.; Marx, D. The shapes of protons in hydrogen bonds depend on the bond length. Chem. Phys. Chem. 2005, 6, 1738–1741. (66) Miller, W. H.; Handy, N. C.; Adams, J. E. Reaction path Hamiltonian for polyatomic molecules. J. Chem. Phys. 1980, 72, 99–112.

24


Page 24 of 39

Page 25 of 39 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 1: Characteristic distances (in ˚ A) and angles (in degree) in the most stable structure of a 2,4-dichlorophenol-trimethylamine complex optimized with various methods. 3ob, NHmod, and NHopt represent respectively SCC-DFTB3 calculations using the 3ob parameter set, using 3ob with a modified N-H repulsive potential, and using 3ob with the N-H repulsive potential optimized in this work. B3LYP and MP2 denote the levels of the calculations with the 6-31G(d,p) basis set. r = |rN − rH | − |rO − rH | and R = |rN − rO |. 3ob NHmod NHopt B3LYP MP2 r 1.119 1.111 0.824 0.800 0.754 R 3.000 2.992 2.740 2.728 2.677 ∠OHN 152.4 152.4 152.1 150.7 151.4

25



Figure 1: (a) 2,4-dichlorophenol-trimethylamine complex handled in this work. (b) Protontrimethylamine complex used to fit the N-H pair repulsive potential.

26


Page 26 of 39

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


Figure 2: Potential energy profiles for the proton transfer reaction along the reaction coordinate r = |rN − rH | − |rO − rH | in gas phase. The minimum at around r = 1.0 is set to zero in each profile. MP2 and B3LYP results almost overlap each other. Results from the SCC-DFTB3 calculations with three different N-H pair repulsive potentials are plotted: DFTB3(3ob) used the original one; DFTB3(3ob + NHmod) used a modified one developed previously; DFTB3(3ob + NHopt) used the N-H repulsive potential refined for the specific system studied in this work. The grid spacing is 0.01 ˚ A.

27



Figure 3: N-H pair repulsive potentials used in the SCC-DFTB3 calculations. [3ob] is the original repulsive potential; [NHmod] is a modified one previously developed; [NHopt] is that constructed in this work.

28


Page 28 of 39

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


Figure 4: The free energy surfaces for the proton transfer reaction along r = |rN −rH |−|rO − rH | and R = |rN − rO | calculated from (a) the SCC-DFTB-MD calculations and from (b) the SCC-DFTB-PIMD calculations. The contour lines are plotted in 2 kcal/mol step. The minimum free energy paths obtained from the corresponding one-dimensional free energy profile calculations along r as the reaction coordinate are plotted with broken lines. (c) The difference between the classical and quantum free energy surfaces. 29



Figure 5: The free energy profile of the proton transfer reaction with R = 2.8 ˚ A obtained from the SCC-DFTB-MD calculations (broken line) and from the SCC-DFTB-PIMD calculations (solid line).

30


Page 30 of 39

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


Figure 6: (a) The free energy profile along the minimum free energy path of the proton transfer reaction calculated from the SCC-DFTB-MD (broken line) and from the SCCDFTB-PIMD (solid line). The minimum of this free energy profile is set to zero for each profile. (b) The difference between the free energy profiles shown in (a). The grid spacing is 0.02 ˚ A.

31



Figure 7: (a) The distance between each bead and the ring polymer centroid G. (b) The force vector acting on a bead (Fk ) and the vector from the bead to the ring polymer centroid (rkG ). The projected force acting on the bead used in the present analysis is obtained by projecting Fk onto rkG /|rkG |.

32


Page 32 of 39

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


Figure 8: Degree of expansion ∆|rH | (see the main text) of the transferring proton ring polymer along the paths with R = |rN − rO | = 2.3, 2.7, and 3.1 ˚ A.

33



Figure 9: Degree of expansion ∆|rH | of the transferring proton ring polymer along the minimum free energy path (broken line), and the projected force acting on the proton ring polymer beads averaged over all the beads and trajectories (solid line). The grid spacing is 0.02 ˚ A.

34


Page 34 of 39

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


Figure 10: Free energy profiles for the proton transfer reaction in a 2,4-dichlorophenoltrimethylamine complex in gas phase along the reaction coordinate r = |rN − rH | − |rO − rH | calculated from the SCC-DFTB-PIMD 10 ps trajectories with various number of polymer beads.

35



Figure 11: Free energy profiles for the proton transfer reaction in a 2,4-dichlorophenoltrimethylamine complex in gas phase along the reaction coordinate r = |rN − rH | − |rO − rH | calculated from the SCC-DFTB-PIMD trajectories of 10 ∼ 50 ps lengths.

36


Page 36 of 39

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


Figure 12: The free energy profiles of the proton transfer reaction at R = 2.8 ˚ A with various timesteps calculated with the SCC-DFTB-PIMD method.

37



Figure 13: Dependence of the barrier height (free energy difference between at r = −0.2 ˚ A and r = 0.9 ˚ A along the free energy profile with R = 2.8 ˚ A) on the timestep used in the SCC-DFTB-PIMD calculations.

38


Page 38 of 39

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


Figure 14: TOC and Abstract graphics.

39


The SCC-DFTB-PIMD method to evaluate a multidimensional

Recommend Documents