Theoretical Study on Excited States of Bacteriochlorophyll a in

Subscriber access provided by UNIV OF TASMANIA

Article

Theoretical Study on Excited States of Bacteriochlorophyll a in Solutions with Density Functional Assessment Masahiro Higashi, Takahiro Kosugi, Shigehiko Hayashi, and Shinji Saito J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp507259g • Publication Date (Web): 25 Aug 2014 Downloaded from http://pubs.acs.org on September 1, 2014

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

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

Page 1 of 48

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

The Journal of Physical Chemistry

Prepared for J. Phys. Chem. B Revised on August 21, 2014

Theoretical Study on Excited States of Bacteriochlorophyll a in Solutions with Density Functional Assessment Masahiro Higashi†*, Takahiro Kosugi§‡, Shigehiko Hayashi§, and Shinji Saito║#* †

Department of Chemistry, Biology and Marine Science, University of the Ryukyus, 1

Senbaru, Nishihara, Okinawa 903-0213, Japan §

Department of Chemistry, Graduate School of Science, Kyoto University, Kyoto 606-8502,

Japan ║

Department of Theoretical and Computational Molecular Science, Institute for Molecular

Science, 38 Nishigo-Naka, Myodaiji, Okazaki, Aichi 444-8585, Japan #

The Graduate University for Advanced Studies (SOKENDAI), 38 Nishigo-Naka, Myodaiji,

Okazaki, Aichi 444-8585, Japan

‡

Present address: Department of Biochemistry, University of Washington, Seattle,

Washington 98105, United States

ACS Paragon Plus Environment

1


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 2 of 48

ABSTRACT: The excited-state properties of bacteriochlorophyll (BChl) a in triethylamine, 1-propanol, and methanol are investigated with the time-dependent density functional theory by using the QM/MM reweighting free energy SCF method. It is found that no prevalent density functionals can reproduce the experimental excited-state properties, i.e. the absorption and reorganization energies, of BChl a in the solutions. The parameter  in the range-separated hybrid functional is therefore optimized to reproduce the differences of the absorption energies in the solutions. We examine the origin of the difference of the absorption energies in the solutions, and find that sensitive balance between contributions of structural changes and solute-solvent interactions determines the differences. The accurate description of the excitation with the density functional with the adjusted parameter is therefore essential to the understanding of the excited-state properties of BChl a in proteins and also the mechanism of the photosynthetic systems.

KEYWORDS: BChl a, QM/MM-RWFE-SCF, TDDFT, Absorption energy, Reorganization energy, Solvent effect


2

Page 3 of 48

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


1. INTRODUCTION Bacteriochlorophyll (BChl) a (Fig. 1) plays important roles in efficient energy transfer of photosynthetic systems as energy carriers in light-harvesting antennas such as FennaMatthews-Olson (FMO) complex, light-harvesting complex 1 (LH1) and LH2.1-8 Thus the excited-state properties of BChl a not only in proteins9-19 but also in solutions,9,19-33 lowtemperature glasses,9,25,31,34 and polymer films10-11 have been experimentally studied for a long time.

Figure 1. Molecular structure of bacteriochlorophyll a. The orientation of x- and y-axis is also shown.

The transition energies to the first excited state, called Qy state from the direction of transition dipole moment (Fig. 1), are widely distributed (up to ~1000 cm−1) in proteins, e.g. FMO complex35 and LH2 antenna.4 In contrast, the transition energies to the Qy state in solutions are rather insensitive to solvent and only slightly sensitive (~100 cm−1) to aliphatic alcohols.24,29,33 For example, the absorption energy to the Qy state of BChl a in 1-propanol is slightly smaller compared with those in methanol and triethylamine. Other excited-state


3


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 48

properties of Qy state have also been investigated. The differences between dipole moments in the ground and Qy states of BChl a are much different in polymer films and proteins.10-14,1718

The reorganization energy for Qy absorption of five-coordinated BChl a in triethylamine is

50% larger than that for Qy emission, while that for Qy absorption of six-coordinated BChl a in a low-temperature triethylamine glass is almost the same as that for Qy emission.31 On the other hand, the transition energies to the second excited state, called Qx state, in solutions strongly depend on the coordination number of the ligand(s), i.e. the solvent molecule(s) coordinated to Mg2+ ion at the center of BChl a.20-21,28-29,31 For example, the absorption energy to the Qx state of five-coordinated BChl a in triethylamine is large compared with those of six-coordinated BChl a in methanol and 1-propanol. However, BChl a in proteins is mostly five-coordinated, and thus the transition energies to the Qx state in proteins are less sensitive than those in solutions. These results suggest that the excited states of BChl a significantly depend on solution and protein environment. In addition to the studies on these fundamental excitedstate properties of BChl a, the excitation energy transfer (EET) dynamics in the lightharvesting antennas have been intensively studied experimentally6,36-39 since the finding of long-lived quantum coherence in the FMO complex with two-dimensional electronic spectroscopy40. Various theoretical studies using model Hamiltonian have also been performed to elucidate the long-lived coherent dynamics.41-48 For example, Ishizaki et al. showed the importance of the finite timescale effects of protein-induced fluctuationdissipation, i.e. non-Markovian interplay between electronic excitation and the surrounding protein, on the EET dynamics.41-42 Recently, several theoretical studies using all-atom molecular dynamics (MD) simulation or normal mode analysis combined with electronic structure calculations have been carried out to understand the EET dynamics in the FMO complex at the atomistic level.49-57 In all these studies, the excitation energies of each BChl a


4

Page 5 of 48

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


are calculated along the MD classical simulations or the normal modes. The obtained EET dynamics are, however, significantly different among these studies because of the difference in the site energies and their fluctuations used in the calculations. As pointed out in Ref. 52, one of the main reasons for the difference between EET dynamics is considered to be due to the difference of the calculated electronic structure, in particular, excited-state properties, of BChl a. Therefore, a proper description of the excited-state properties, such as the excitation energies and reorganization energies, of BChl a is essential to the understanding of the EET dynamics. Here, we focus on the excited-state properties of BChl a in solutions. Although there are many theoretical studies on the electronic structures of the excited states of BChl a,31,35,5877

the excited-state properties in solutions, such as the solvent (in)sensitivty to the Qy and Qx

absorption energies, have not been theoretically investigated, to the best of our knowledge. Three solvents, triethylamine, 1-propanol, and methanol, were used in the present calculations because the experimental data of the reorganization energies of BChl a in triethylamine are available31 and the Qy and Qx states of aliphatic alcohols have been well studied.24,29,33 We examined the excited-state properties of BChl a in the three solutions with quantum mechanical and molecular mechanical (QM/MM) method.78-80 To describe the electronic structure in solutions at ambient temperature, we employed a QM/MM free energy geometry optimization methods called the QM/MM reweighting free energy self-consistent field (QM/MM-RWFE-SCF) method.81-82 The QM/MM-RWFE-SCF method is a combination of two QM/MM free energy geometry optimization methods: a QM/MM free energy optimization method based on a mean field approximation developed by Yamamoto83 and a reweighting update scheme for statistical ensemble of environmental conformation introduced by Hu et al.84 The QM/MM-RWFE-SCF method enables us to optimize the


5


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

Page 6 of 48

geometry of the QM molecule effectively on the free energy surface defined by the surrounding MM environment obtained from MD simulation. For an electronic structure method to describe the QM molecule, we used density functional theory (DFT) and time-dependent density functional theory (TDDFT). Although the TDDFT is a cost-effective powerful tool for calculating the excited states of relatively large molecules (up to ~100 atoms), the accuracy strongly depends on the density functional. We examined various types of density functionals, including pure, global hybrid, and rangeseparated hybrid ones, and found that none of functionals can reproduce the excited-state properties of BChl a in solutions. Thus, as the previous studies,73,85-93 we adjusted the parameter  in the range-separated hybrid functionals to reproduce the relative absorption energies to the Qy state in the three solutions. We analyzed why the conventional functional cannot reproduce the excited-state properties by comparing the result with the adjusted parameter to that with the original one. The effect of hydrogen bonds between alcohol molecules and two conjugated carbonyl groups of BChl a, 3-acetyl and 131-keto groups (Fig. 1) was also analyzed to investigate the solvent (in)sensitivties to the absorption energies. The organization of the article is as follows. In the next section, we describe the theoretical methods employed here. The computational details of the methods are given in section 3. In section 4, we present and discuss the results of the calculations, and the conclusions are summarized in section 5.


6

Page 7 of 48

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


2. THEORETICAL METHODS We calculated the absorption energies and reorganization energies of BChl a in triethylamine, 1-propanol, and methanol by using the QM/MM-RWFE-SCF method.81-82 In the QM/MM-RWFE-SCF calculation, BChl a and ligand(s), i.e. solvent molecule(s) coordinated with Mg2+ ion, were treated quantum mechanically. The geometry of the QM subsystem was optimized in the presence of thermally averaged distribution of solvent molecules described by the MM force field. We briefly describe the QM/MM-RWFE-SCF method here, though the details of the QM/MM-RWFE-SCF method are found elsewhere.81-82 The free energy difference between a target state (Ψ,R ) , where  and R are the electronic wave function and coordinates in the QM region, and a reference thermal equilibrium state (Ψ ref , R ref ) obtained from a MD simulation, is defined as F [,  ref , R, R ref ]  E QM [, R ]  E QM [ ref , R ref ]   1 ln exp     E QM-MM [, R, X ]  E QM-MM [ ref , R ref , X ]

(1)  ref , R ref

where E QM and E QM-MM are the energy of the QM subsystem and the QM-MM interaction energy, respectively, and   1 / kBT is the inverse temperature. The X is the atomic coordinates in the MM region, and

 ref ,R ref

denotes the ensemble average over the thermal

distribution of MM conformations with the QM wave function and the coordinates fixed at (Ψ ref , R ref ) . The QM-MM interaction energy E QM-MM can be divided into electrostatic (ES)

and non-electrostatic (non-ES) parts. The latter part consists of van der Waals and valence interactions and, thus, is generally expressed in terms of force fields, QM-MM QM-MM E QM-MM  , R, X   EES (R, X) . , R, X  Enon-ES

(2)


7


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

Page 8 of 48

We adopt a site-site representation of the QM-MM ES interaction,83,94-95 N QM

E

QM-MM ES

 , R, X    qa    a  R, X 

(3)

a

where N QM is the number of QM atoms. The qa is the partial charge on the QM atomic site a:

qa     qˆa 

(4)

where qˆa is the population operator that generates the partial charge qa , and  a is the ES potential from the MM region, a   A

qAMM Ra  X A

(5)

where qAMM is the effective charge of MM atom A. Note that by adopting this representation, QM-MM EES can be regarded as a function of (q, R,Φ) , which drastically reduces the

computational cost. The average ES potential Φ at (q, R) is obtained as

a

q ,R



 a exp     E QM-MM (q, q ref , R, R ref , X )  exp     E

QM-MM

(q, q ref , R, R ref , X ) 

q ref , R ref

(6)

q ref , R ref

where E QM-MM (q, qref , R, R ref , X)  E QM-MM (q, R, X)  E QM-MM (qref , R ref , X) , and the average wave function in the presence of Φ is obtained with the electronic structure calculation. Note that the electronic structure calculation and Eq. (6) have to be solved in a self-consistent manner because Φ depends on the charges q[] . The QM geometry can be optimized on the free energy surface by using one MD trajectory with the QM charges and coordinates fixed at (q ref , R ref ) because the analytical free energy gradient is available by differentiating Eq. (1) with respect to R. However, Eqs. (1) and (6) are valid only if the difference between (q ref , R ref ) and (q, R ) is small because of


8

Page 9 of 48

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


the finite length of the MD trajectory. Thus, it is necessary to perform the sequential sampling,84 where MD simulations are iteratively repeated to update the ensemble for the renewed (q ref , R ref ) , which is taken from the optimized geometry in the preceding step, until the difference between the reference (q ref , R ref ) and optimized variables (q, R ) is sufficiently small. The Q ( = x, y) absorption energy of BChl a in solutions was calculated at the geometry optimized in the S0 ground state,

Eabs (Q )  E EEQM (Q ; R0 , Φ 0 )  E EEQM (S0 ; R0 , Φ 0 )

(7)

QM-MM where E EEQM is the electrostatically embedded QM energy,95 E EEQM  E QM  EES . The

R0 and Φ

0

are the geometry and ES potential optimized in the ground state. The

reorganization (free) energy for Qy absorption, abs, was calculated as

abs  Eabs (Qy )  F (S0 ; R0 )  F (Qy ; R y )

(8)

where F (S0 ; R 0 ) and F (Qy ; R y ) are the free energies of the ground and Qy states at the geometry optimized on each free energy surface, and the difference of the second and third terms in Eq. (8) was evaluated by using the QM/MM free energy perturbation method. The reorganization energy for Qy emission, emi, was similarly obtained as

emi   Eemi (Qy )  F (Qy ; R y )  F (S0 ; R0 )

(9)

Eemi (Qy )  E EEQM (Qy ; R y , Φ y )  E EEQM (S0 ; R y , Φ y ) .

(10)

where

The importance of inclusion of solvent molecules around a solute molecule into the QM region for the excited-state properties has been reported.74,96-97 It is considered that the inclusion of alcohol molecules forming hydrogen bonds with two conjugate carbonyl groups of BChl a, 3-acetyl and 131-keto groups (Fig. 1) affects the excited-state properties, though the solvent molecules are described by the MM force field in the above procedure. Therefore,


9


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

Page 10 of 48

we added the quantum mechanical correction to the absorption energy calculated with the QM/MM-RWFE-SCF method in alcohol solutions by using the subset of the ensemble in the QM/MM-RWFE-SCF procedure, wQC QC Eabs (Q )  Eabs (Q )  Eabs (Q )

(11)

where QC large small Eabs (Q )  Eabs (Q ; R0 , X)  Eabs (Q ; R0 , X)

subset

.

(12)

large Eabs (Q ; R0 , X) is the Qabsorption energy calculated with the large QM subsystem with

alcohol molecules making hydrogen bonds with the two conjugated carbonyl groups into the small (Q ; R0 , X) is that calculated with the QM subsystem consisting of QM region, while Eabs

BChl a and ligand(s) molecules. There are no hydrogen bonds between the carbonyl groups QC and solvent molecules in triethylamine solution and thus Eabs (Q )  0 .

We used the DFT and TDDFT methods for the QM subsystem. One of the recent great advances in (TD)DFT is the introduction of two range-separated hybrid schemes: longrange-corrected (LC)98-99 and Coulomb-attenuating method (CAM)100 schemes. In these schemes, the two-electron repulsion operator 1/r12 in the exchange functional is divided into short- and long-range parts using the error function,

1 1     erf(  r12 )    erf(  r12 )   , r12 r12 r12

(13)

with three parameters,  , and . The DFT exchange functional is included through the first term, and the long-range orbital-orbital exchange interaction is described using the HartreeFock exchange integral. The distance-dependent ratio of these two terms is determined by the range-separated parameter . Based on Eq. (13), various LC and CAM functionals have been developed, e.g., LC-BOP,99,101 LC-BLYP,101 and CAM-B3LYP,100 where the LC is a special case ( = 0 and  = 1) of the CAM. From the definition, as  → 0, LC-BLYP and CAM-


10

Page 11 of 48

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


B3LYP approach BLYP and B3LYP functional. (Strictly speaking, CAM-B3LYP with = 0 is different from B3LYP due to the factor for Becke’s 1988 gradient correction in B3LYP.) The parameters in Eq. (13),  , and ,have been determined to minimize the mean error of a benchmark database, which generally consists of experimental data or high-level calculation results of atoms and small molecules (up to ~10 atoms). For example, Hirao and co-workers first determined  = 0.33 in the LC scheme to minimize the mean absolute error of the optimized bond lengths of homonuclear diatomic molecules for the first- to third-row atoms except rare-gas atoms.101 Then, they obtained  = 0.47 to minimize the root mean square percent error for the atomization energies of the G2 set102 (148 molecules).103 The parameters of CAM-B3LYP,  = 0.19 and  = 0.46, were determined to yield the least error for the atomization energies of the small G2 set (53 molecules) with  = 0.33.100 Recent studies show, however, that such parameters are not appropriate for describing excited states of large molecules.73,85-93 For instance, it was reported89 that the vertical transition energies to the valence excited states of large molecules are well reproduced with the LC-PBE functional104-105 with  = 0.20, which is smaller than the originally optimized value,  = 0.40. Kityk calculated the absorption and fluorescence spectra of two organic heterocyclic dyes in several solutions with LC-BLYP in the combination with the polarizable continuum model (PCM) and found that  = 0.14-0.15 and  = 0.21-0.23 are the best estimates for absorption and fluorescence energies, respectively.93 Smaller values of are recommended for large molecules compared with the originally optimized ones, although the best parameter  depends on the systems and properties to be optimized. As shown below, we optimized the parameter  to reproduce the relative absorption energy to the Qy state of BChl a in solution properly.


11


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

Page 12 of 48

3. COMPUTATIONAL DETAILS The QM subsystems in the QM/MM-RWFE-SCF calculations consist of BChl a and ligand(s) coordinated to Mg2+ ion. It is considered from the experimental Qx absorption energies in solutions that the Mg2+ ion of BChl a in triethylamine is five-coordinated while those in 1-propanol and methanol are mostly six-coordinated.20-21,28-29,31 According to these results, we added one and two ligands to the QM subsystem for triethylamine solution and for 1-propanol and methanol solutions, respectively. To reduce computational cost, we replaced long alkyl chains with methyl groups; two long alkyl chains of BChl a, 8-ethyl and 17propionate groups (the latter is usually esterified with phytol), were replaced with methyl ones. However, the effect of the phytyl tail on the coordination number of Mg, which affects the Qy absorption energy slightly, has been experimentally known.29 For ligand(s), triethylamine was replaced with trimethylamine, and 1-propanol was replaced with methanol. As a result, the QM subsystem in 1-propanol is the same as that in methanol (Fig. 2).

Figure 2. Structures of QM subsystems, BChl a with ligand(s): in triethylamine (left) and in 1-propanol and methanol (right).

The QM/MM systems for triethylamine (TEA), 1-propanol (PRO), and methanol (MET) solutions consist of 1 solute and 706 triethylamine, 1207 1-propanol, and 2001 methanol molecules, respectively. The box lengths of the cubic boxes are 54.7, 53,5, and 51.8 Å, respectively, and the periodic boundary conditions were used. The OPLS-AA force field modified by Kahn and Bruce106-107 was used for ligand and solvent molecules. The van der


12

Page 13 of 48

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


Waals interactions between solute and solvent molecules were also treated with the modified OPLS-AA force field. Although the results of the triethylamine and methanol solutions calculated with the modified OPLS-AA force field are similar to those with the original OPLS-AA one, the description of the Qy absorption energy for the 1-propanol solution by the modified OPLS-AA force field is better than that with the original OPLS-AA one. We first carried out the MD simulations of triethylamine and methanol solution systems with MM force field for all the molecules to obtain the initial geometries of the QM subsystem for the QM/MM-RWFE-SCF calculations. For triethylamine solution, we added the trimethylamine ligand in the −z region (Fig. 2). Notably, our calculations on the absorption energies of BChl a with the ligand in the gas phase indeed show that the side of ligation little affects the absorption energies compared with the solvent effects, e.g. ~30 and ~20 cm−1 for the Qy absorption energies with B3LYP and CAM-B3LYP ( = 0.20). The AMBER GAFF force field108 was used for the intramolecular potential energy function and partial charges of BChl a molecule. In addition, weak harmonic restraints were imposed on the distance between Mg and ligating atoms to avoid ligand exchanges and to keep the fiveand six-coordinations for Mg in triethylamine and methanol solutions. Long-range ES interactions were treated with the particle mesh Ewald method.109 Bonds involving hydrogen were constrained by using the SHAKE method.110 Equilibrium MD simulations at 300 K were carried out for 500 ps. Then the conventional QM/MM geometry optimizations94 were carried out for the last snapshots of MD simulations. In the QM/MM geometry optimizations, atoms within 20 Å from the Mg were allowed to move, while the outer atoms were fixed. The long-range ES interaction in the QM/MM calculation was smoothly truncated at 15 Å by using the tapering method implemented in the TINKER program.111 The optimized QM geometries and partial charges were used as the initial reference of the QM/MM-RWFE-SCF


13


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

Page 14 of 48

calculations. The geometry and partial charges optimized in methanol solution were employed as the initial reference for 1-propanol solution system. In the QM/MM-RWFE-SCF calculations, the sampling MD simulations were carried out with the QM geometries fixed by using the SHAKE method. A trajectory for 6 ns was carried out at each iteration step of the sequential sampling for the triethylamine and methanol solution systems. A longer trajectory with 12 ns was performed in each interaction step for 1-propanol solution because the hydrogen bond dynamics is slow.112 From these trajectories, 30,000 configurations were taken as the MM ensemble; i.e., the MM conformational samples were taken at every 300 fs for triethylamine and methanol and 600 fs for 1-propanol, respectively. The convergence criterion of the QM geometry optimization at each step of the sequential sampling was set to be 5.0 × 10−4 Hartree/Bohr for the largest component of the QM gradient. The sequential samplings were repeated until all the rootmean-square deviations of QM optimized geometries from those at the preceding steps are less than 1.0 × 10−3 Å in three successive QM/MM-RWFE-SCF optimization steps. Most of the QM geometry optimizations were converged within 10 sequential samplings. The last MD trajectories for 6 ns and 12 ns for triethylamine and methanol, and 1-propanol solutions were employed for further analysis. The quantum mechanical correction to the absorption energy, Eq. (12), was calculated by averaging the differences between the absorption energies with the large and small QM subsystems at 100 configurations, which were taken from 30,000 configurations used in the last QM/MM-RWFE-SCF optimization step. We defined that a hydrogen bond between the carbonyl group and an alcohol molecule is formed when the distance ROHO < 2.9 Å, ROOH < 3.5 Å, and the angle OOH HO < 30 degree,113-114 where the distances are taken from the first minimum of the radial distribution functions.


14

Page 15 of 48

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


We used the DFT and TDDFT methods with several types of density functionals in QM/MM calculations: two pure functionals, BLYP and PBE, two global hybrid functionals, B3LYP and PBE0, and two range-separated hybrid functionals, CAM-B3LYP and LC-BLYP. As mentioned above, two parameter values,  = 0.33 and 0.47, were proposed for LCBLYP.101,103 Our preliminary calculations on the absorption energies in three solutions showed that the results with  = 0.33 are better than those with  = 0.47. (In fact, as shown below, a smaller value is preferred.) Therefore, only the results with  = 0.33 are shown hereafter. We also tested the HF/CIS method for comparison. We employed the 6-31G(d,p) basis set. We also tested the 6-311G(d,p) and 6-31+G(d,p) basis sets, but no significant improvement was found, in fact, worse results were found in some cases, e.g. reorganization energies. We adopted the RESP scheme94,115 to describe partial charge qa . The absorption energies and reorganization energies of BChl a in the gas phase were calculated without ligand and solvent molecules. The reorganization energies in the gas phase were directly obtained as the difference between the energies at the optimized geometries in the ground and excited states.31 All the MD simulations were carried out using the AMBER 11 program package.116 The electronic structure calculations were performed using the GAMESS program package,117 where we implemented the QM/MM potential energy and QM/MM-RWFE-SCF routines.


15


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

Page 16 of 48

4. RESULTS AND DISCUSSION Tables 1 and 2 summarize the calculated absorption energies to the Qy and Qx states with the quantum mechanical correction for hydrogen bond interactions with the experimental results. The experimental Qy absorption energies in triethylamine (dielectric constant = 2.4), 1-propanol ( = 20.5), and methanol ( = 32.7) solutions are 12991, 12870, and 12974 cm−1, respectively;29 The Qy absorption energy in 1-propanol is slightly lower than those in methanol and triethylamine, Eabs(Qy; PRO) < Eabs(Qy; MET) ~ Eabs(Qy; TEA). The Table 1. Calculated Qy Absorption Energies (in cm−1) of BChl a in the Gas Phase, Triethylamine, 1-Propanol, and Methanol. Values in Parentheses are the Differences from the Absorption Energies in Triethylamine. Method

Gas

Triethylamine

BLYP

14590

14580

13940a

(−640)

13820a

(−760)

PBE

14750

14730

14190a

(−540)

14000a

(−730)

B3LYP

15180

15180

14890

(−290)

14820

(−360)

PBE0

15340

15350

15090

(−260)

15040

(−310)

CAM-B3LYP ( = 0.33)

14370

15140

15470

(330)

15820

(680)

LC-BLYP ( = 0.33)

16880

17980

19370b

(1390)

19690b

(1710)

HF/CIS

29220b

30260b

31730b

(1470)

32030b

(1770)

CAM-B3LYP ( = 0.20)

14580

14680

14570

(−110)

14650

(−30)

LC-BLYP ( = 0.18)

14070

14050

14030

(−20)

14050

(0)

12991

12870

(−121)

12974

(−17)

Experimentc

1-Propanol

a

The first excited state composed of mixtures of Qy and Qx states.

b

The second excited state.

c

Ref. 29.

Methanol


16

Page 17 of 48

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


experimental Qx absorption energies in triethylamine, 1-propanol, and methanol solutions are 17323, 16454, and 16461 cm−1, respectively;29 i.e., Eabs(Qx; TEA) > Eabs(Qx; PRO) ~ Eabs(Qx; MET). We assigned to the Qy and Qx states based on the direction of the transition dipole moment. In most cases, as the experimental results, the first and second excited states were easily identified as the Qy and Qx states, respectively. The dominant contribution to the Qy transition is the HOMO → LUMO transition, while that to the Qx transition is the HOMO−1 Table 2. Calculated Qx Absorption Energies (in cm−1) of BChl a in the Gas Phase, Triethylamine, 1-Propanol, and Methanol. Values in Parentheses are the Differences from Absorption Energies in Triethylamine. Method

Gas

Triethylamine

BLYP

16410a

15660

14360b

(−1300)

14300b

(−1360)

PBE

16540a

15780

14540b

(−1240)

14440b

(−1340)

B3LYP

18230

17580

16160

(−1420)

16070

(−1510)

PBE0

18700

18070

16730

(−1340)

16660

(−1410)

CAM-B3LYP ( = 0.33)

18930

18650

17570

(−1080)

17670

(−980)

LC-BLYP ( = 0.33)

19670

19220

18280c

(−940)

18350c

(−870)

HF/CIS

26950c

26620c

26060c

(−560)

26100c

(−520)

CAM-B3LYP ( = 0.20)

18510

17900

16670

(−1230)

16670

(−1230)

LC-BLYP ( = 0.18)

17700

16900

15780

(−1120)

15780

(−1120)

17323

16454

(−869)

16461

(−862)

Experimentd

1-Propanol

a

The third excited state.

b

The second excited state composed of mixture of Qy and Qx states.

c

The first excited state.

d

Ref. 29.

Methanol


17


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 18 of 48

→ LUMO transition. With the pure density functionals, i.e., both BLYP and PBE, however, the Qx state is the third excited state in the gas phase because a dark excited state represented by the HOMO−2 → LUMO transition lies between the Qy and Qx states. In addition, the assignments of the Qy and Qx states are difficult in alcohol solutions with the pure density functionals because the excitation energy differences between the first and second excited states are small (at most ~500 cm−1) and these excitations are characterized as the mixture of the HOMO → LUMO and HOMO−1 → LUMO transitions. The mixing of the HOMO → LUMO and HOMO−1 → LUMO transitions is also found in the previous studies on the excitation energies of BChl a and other chlorophylls calculated with pure density functionals.77,118 Charge transfer excited states from an alcohol to BChl a are often found near the two states in the present calculation with the large QM subsystem with the pure density functionals, though such states are not found with the other hybrid functionals. We found that the order of the Qy and Qx states in the alcohol solutions is changed in the results with the LC-BLYP ( = 0.33) functional due to the large decrease of the dipole moment in the Qy state (Table 3, see also below). The same change is also found in the HF/CIS calculations. The orders of the Qy and Qx states in the LC-BLYP and HF/CIS calculations are not in agreement with the experimental results. All the calculated Qy absorption energies are overestimated, as the previous studies using the TDDFT and HF/CIS methods.52,61-62,64,71-73,76 Therefore, we focus on the energy differences in three solutions rather than the absorption energies themselves. However, none of the density functionals employed here can reproduce the order of Qy absorption energies in three solutions; Eabs(Qy; MET) < Eabs(Qy; PRO) < Eabs(Qy; TEA) in the pure and global hybrid density functionals, whereas Eabs(Qy; TEA) < Eabs(Qy; PRO) < Eabs(Qy; MET) in the rangeseparated hybrid functionals and the HF/CIS method.


18

Page 19 of 48

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


Unlike the Qy excitation energy difference, the Qx excitation energy difference between triethylamine and alcohol solutions is qualitatively reproduced in all the calculations. However, the difference between the excitation energies in the two alcohol solutions is not well reproduced; Eabs(Qx; MET) < Eabs(Qx; PRO) in the pure and global hybrid functionals and

Eabs(Qx; PRO) < Eabs(Qx; MET) in the range-separated hybrid functionals and the

HF/CIS method.

Table 3. Differences between Dipole Moments (in Debye) of BChl a in the Ground and Qy States, d y , in the Gas Phase, Triethylamine, 1-Propanol, and Methanol at the Ground-State Optimized Geometries. Method

Gas

Triethylamine

1-Propanol

Methanol

BLYP

0.30

0.41

0.48

0.52

PBE

0.31

0.41

0.49

0.50

B3LYP

0.31

0.46

0.60

0.65

PBE0

0.33

0.54

0.66

0.74

CAM-B3LYP ( = 0.33)

1.99

3.70

4.54

5.01

LC-BLYP ( = 0.33)

7.95

7.91

8.27

8.38

HF/CIS

9.97

9.86

9.16

8.92

CAM-B3LYP ( = 0.20)

0.68

1.27

1.72

2.09

LC-BLYP ( = 0.18)

0.65

1.03

1.60

1.79

The above results demonstrate that the calculated absorption energies in solutions significantly depend on the density functional. Such dependencies of the density functional on the absorption energies of BChl a and other chlorophylls are also found in the previous studies.77,118-119 Next, we examined the dipole moments in the ground, Qy, and Qx states. The


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

Page 20 of 48

differences between dipole moments in the ground and Q states, d  d  d0 , at the ground-state optimized geometry are shown in Table 3 and Table S1 in the Supporting Information. The dipole moment direction of BChl a in the ground state is almost parallel to y-axis due to the presence of polar 3-acetyl group in the +y region and 131-keto and 132methylcarboxy groups in the −y region (Fig. 1). The dipole moments in the Qy and Qx states are smaller than that in the ground state. The difference between the dipole moments in the ground and Qy states, d y , strongly depends on the density functional. The differences with the pure and global hybrid functionals are rather small (

Theoretical Study on Excited States of Bacteriochlorophyll a in

Recommend Documents