Subscriber access provided by Stony Brook University | University Libraries
Spectroscopy and Excited States
A First-Order Interacting Space Approach to Excited-State Molecular Interaction: Solvatochromic Shift of p-Coumaric Acid and Retinal Schiff Base Kazuma Yanai, Kazuya Ishimura, Akira Nakayama, and Jun-ya Hasegawa J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01089 • Publication Date (Web): 11 Jun 2018 Downloaded from http://pubs.acs.org on June 20, 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 45 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
A First-Order Interacting Space Approach to Excited-State Molecular Interaction: Solvatochromic Shift of p-Coumaric Acid and Retinal Schiff Base Kazuma Yanai,† Kazuya Ishimura,‡,* Akira Nakayama,† and Jun-ya Hasegawa†,* †
Institute for Catalysis, Hokkaido University, Kita 21, Nishi 10, Kita-ku, Sapporo 001-0021,
Japan ‡
Department of Theoretical and Computational Molecular Science, Institute for Molecular
Science, 38 Nishigo-Naka, Myodaiji, Okazaki 444-8585, Japan * Corresponding author, e-mail:
[email protected],
[email protected]
KEYWORDS Excited-state molecular interactions, First-order interacting space, Energy decomposition analysis, Solvatochromism
ABBREVIATIONS AS, active site; bR, bacteriorhodopsin; CASPT2, complete active space second-order perturbation theory; CIS, configuration interaction singles; CT, charge transfer; DPSB,
ACS Paragon Plus Environment
1
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 2 of 45
deprotonated retinal Schiff base; FOIS, first-order interacting space; LMO, localized MO; MD, molecular dynamics; p-CA, para-coumaric acid; [p-CAp]-, phenolate form of p-CA anion; [pCAc]-, carboxylate form of the p-CA anion; PSB, protonated retinal Schiff base, SAC-CI, symmetry adapted cluster-configuration interaction; sQM, secondary QM;.
ABSTRACT
A triple-layer QM/sQM/MM method was developed for accurately describing the excited-state molecular interactions between chromophore and the molecular environment (J. Hasegawa, K. Yanai, and K. Ishimura, ChemPhysChem 2015, 16, 305). A first-order-interaction space (FOIS) was defined for the interactions between QM and secondary QM (sQM) regions. Moreover, configuration interaction singles (CIS) and its second-order perturbation theory (PT2) calculations were performed within this space. In this study, numerical implementation of this FOISPT2 method significantly reduced the computing time, which realized application to solvatochromic systems, para-coumaric acid in neutral (p-CA) and anionic forms in aqueous solution, retinal Schiff base in methanol (MeOH) solution, and bacteriorhodopsin (bR). The results were consistent with the experimentally observed absorption spectra of the applied systems. The QM/sQM/MM result for the opsin shift was in better agreement to the experimental result than that of the ordinary QM/MM. A decomposition analysis was performed for the excited-state molecular interactions. Among the electronic interactions, charge-transfer (CT) effect, excitonic interaction, and dispersion interaction showed significant large contributions,
ACS Paragon Plus Environment
2
Page 3 of 45 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
while electronic polarization effect presented only minor contribution. Furthermore, the result was analyzed to determine the contributions from each environmental molecule and was interpreted based on the distance of the molecules from the π-system in the chromophores.
ACS Paragon Plus Environment
3
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 45
1. Introduction In biological systems, interactions between light and molecules are highly relevant to their functions. The wavelength of the interacting light is crucially important in color vision1-4 and biological emission5-6 to detect/illuminate particular color of light for their purposes. For example, our color recognition is based on three visual pigments with different photo-absorption wavelength.1-4 Similarly, artificial molecular systems are required to tune absorption/emission wavelength in the fields of photochemistry7-8 and materials science.9 To understand and rationally control the absorption/emission energy, it is necessary to calculate excitation energy in a complex molecular environment. There are several approaches to treat molecular interactions for calculating absorption/emission spectrum (for review, see references
10-13
) Hybrid approaches with different combination of methods were developed,
including dielectric models,11, 14 RISM-SCF-SEDD theory,15-16 MMpol models,17-18, embedded fragment potential,19 frozen density embedding,20 and effective fragment potential.21 We have also focused on electronic-structure calculations to obtain the photo-absorption/emission energy of chromophore/emitter (solute) in molecular environments (solvent) such as proteins and solutions. Correlated wave function and symmetry-adapted cluster-configuration interaction (SAC-CI) methods22-24 were combined with molecular mechanics methods25 in a hybrid quantum mechanics/molecular mechanics (QM/MM) framework.26 With this approach, spectral tuning mechanisms were studied for several biological systems.25, 27-29 The QM/MM-type calculations were extensively applied to analyze the mechanism of spectral tuning in retinal proteins and successfully interpreted specific roles of molecular environment.25,
27, 30-45
On the other hand, our QM/MM calculation was applied to the
ACS Paragon Plus Environment
4
Page 5 of 45 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
solvatochromic shift in deprotonated retinal Schiff base (DPSB) but failed in reproducing the experimentally observed solvatochromic shift.46 A simple fixed point-charge model failed to explain bathochromic shift when DPSB are introduced to bacteriorhodopsin (bR) environment. This suggests that the interaction between the solute and solvent should be described in a more accurate and reliable way. The case for rhodopsin has been addressed for many years and improvement was achieved, such as incorporating polarization contributions into INDO/S,30-31 CASSCF/CASPT239 and TDDFT45 calculations. The different multipole expansion approaches were adopted, and the results were also different; some showed bathochromic shifts (-0.34 eV,3031
about -0.14 eV,45), and the LoProp result depended on the systems (-0.06 eV to +0.03 eV39).
A triple-layer QM/QM/MM approach was also developed, and DFT calculation was used for the secondary QM region for describing electronic polarization effect.41 The contribution of the polarization effect was much smaller (0.03 eV). Similarly, we extended a QM/MM method to better describe the orbital delocalization and excitonic interactions in excited states.46-48 In this approach, a secondary QM (sQM) region was introduced to improve the interactions between solute and nearby solvent molecules. A (QM+sQM)/MM calculation was performed and an ONIOM49-type correction to the excitation energy was accomplished. Hereinafter, this calculation scheme is named “QM/sQM/MM”. For the calculation of the QM+sQM region, configuration interaction singles (CIS) was adopted, and the orbital delocalization and excitonic interactions were considered at the simplest level. This correction worked well; in the DPSB case, excitation energy in bR relative to that in MeOH solution was improved by 0.11-0.16 eV.4647
In fluorescent proteins, fluorescence energy of mKO relative to that of GFP was improved by
0.05 eV.29
ACS Paragon Plus Environment
5
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 6 of 45
Further improvement on the description of the interactions between QM and sQM regions was developed; with configuration interactions singles (CIS) as a reference wave function, the first-order interacting space was solved using the second-order perturbation theory (FOISPT2). Our idea is to include excited-state molecular interactions in a consistent way in terms of the perturbation-theoretical point of view. A theoretical formalism and preliminary application were already reported on a previous communication48 (a brief summary is given in the next section). However, the inability to calculate realistic systems was encountered. Consequently, applications were performed only for the solvatochromic shift of s-trans-acrolein with 13 H2O molecules and methylenecyclopentene in 22 H2O molecules. In this study, this problem has been resolved by extensive revisions in the source code that take QM+sQM region with more than 1000 atoms. A brief description of the theoretical model (see ref [48] for original manuscript) and computational details were given in sections 2 and 3, respectively. The most effective acceleration was achieved by improved performance in the integral transformation code, which was explained in subsection 2.2. In section 4, the result of the application to solvatochromic shift in para-coumaric acid in neutral (p-CA) and anion forms was presented. In section 5, the results of the spectral tuning mechanism of retinal Schiff base was discussed. In these applications, molecular interactions in excited states were decomposed into contributions from each solvent. A conclusion was given in the final section.
2. Methods 2.1 A triple-layer hybrid QM/sQM/MM method with first-order interacting space (FOIS) approach to excited-state molecular interaction
ACS Paragon Plus Environment
6
Page 7 of 45 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
Our aim is to accurately calculate the excited states of solute molecule surrounded by solvent. We specify “solvent” to include molecular environment such as proteins. Ordinary QM/MM calculations consider the electrostatic interaction between the QM and MM for the excited-state calculation of the QM part. The electrostatic interaction is a Coulombic interaction between the electron density of solute and electrostatic potential from the MM charges. However, as shown in previous studies,14-15, 17, 18, Söderhjelm, 2009 #91, 21, 29-31, 41, 46-48, 50 significant contributions originate from electronic interactions between solute and solvents. At the secondorder perturbation theoretical analysis,48 excitonic interactions, orbital delocalization, electronic polarization, and dispersion interactions appear in the perturbation correction for the excitation energy of singly excited states. These interactions were effectively considered by setting a sQM region between the QM and MM regions; see Figure 1a.29, 41, 46-48 The size of the sQM region is defined by a parameter, rsQM , the atom-atom distance between QM and sQM regions. If a solvent molecule (or an amino acid unit, -CO-CHR-NH-, in proteins) has an atom within rsQM distance from the QM atoms, the solvent molecule is included in the sQM region.
ACS Paragon Plus Environment
7
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 8 of 45
(a) A three-layer model
(b) Fragmentation for LMO Frag. 4
α
α
Frag. 3 Frag. 6
α
Frag. 5 Frag. 1
Frag. 2
Figure 1. (a) Three-layered QM/sQM/MM computational model. (b) Fragmentation of a protein environment for the MO localization. The dotted lines represent the boundaries between each fragment.
In the present QM/sQM/MM scheme, the effect of the sQM region will be introduced by an ONIOM49-type correction. With a certain rsQM value, the energy correction, ∆EY ( rsQM ) , is evaluated as follows:
(
)
(
∆EY ( rsQM ) = EY QM + sQM ( rsQM = d Å ) / MM − EY QM + sQM ( rsQM = 0 Å ) / MM
)
(1)
ACS Paragon Plus Environment
8
Page 9 of 45 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
Based on the equation, rsQM = 0 Å means that there is no sQM region, identical to a standard QM/MM approach. With this correction included, the total energy for I -th state is
(
)
E XI ,Y ( rsQM ) = E XI QM + sQM ( rsQM = 0 Å ) / MM + ∆EYI ( rsQM ) .
(2)
X and Y denote the computational method. The QM+sQM region can be set to be moderately large, e.g., around 1000 atoms for retinal Schiff base with rsQM = 7 Å . Therefore, a possible choice for X and Y would be a highly-accurate method and a fast method with moderate accuracy, respectively. A formula for the excitation energy between I -th and J -th states is readily derived as
(
)
EX XI −,YJ ( rsQM ) = EX XI − J QM + sQM ( rsQM = 0 Å ) / MM + ∆EX YI − J ( rsQM ) ,
(3)
where excitation energies, EX XI −,YJ , EX XI − J , and ∆EX YI − J , are obtained by subtracting each term of eq. 2 for I -th and J -th states. This study employed SAC-CI and CIS + FOISPT2 methods for X and Y, respectively. FOIS is described briefly subsequently. At our starting point, let us assume that, at the starting point, molecular orbitals localized within each solute and solvent and CIS coefficients
{c } 0 ai
∑
ai
for solute were provided by pre-steps. As described in previous studies,48, 51 a commutator,
Hˆ , cai0 Sˆai , is defined, where Hˆ and Sˆai are the Hamiltonian and single excitation (spin
adapted) operators, respectively. Indices a and i denote occupied and unoccupied orbitals, respectively. This commutator includes information about excitation energy.48,
∑
ai
Hˆ , cai0 Sˆai HF
51
When
is projected onto CIS = ∑ HF Sˆai† cai0 , the CIS excitation energy is obtained.
ACS Paragon Plus Environment
9
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 45
We obtain FOIS for the excited-state molecular interactions by projecting out the excitation within a solute:
(1 − Pˆ ) ∑ Hˆ , c 0
ai
frag type
ˆ ˆA ai S ai HF = ∑ ∑ Ωt HF
(4)
A t =1,6
The projection operator Pˆ0 = ∑ I ∈solute I I
is the summation of the excitations within a solute.
Due to the localized molecular orbitals (LMOs), excitation operators that are accompanied by two-electron repulsion integrals such as (φbAφ jB | φa0φi0 ) (MOs b and j belong to different solvent molecules, A ≠ B ) can be neglect. Those integrals do not affect the calculated energy because the LMO product of two different LMOs is very small.52 We used this FOIS {Ωˆ tA HF } to derive the excitation operators for the excited-state molecular interactions. It is composed of six types ( t ) of excitation operators: three kinds of single excitation operators and three kinds of double excitation operators.48 The three single excitations are interpreted as (i) orbital delocalization of solute, (ii) that of solvent, and (iii) excitonic coupling between the solute and solvent. The three double excitations are (iv) electronic polarization of the solute, (v and vi) dispersion interactions in the excited state; see Figure 2 for the schematic illustration of the six kinds of excitation operators.) Owing to the commutator representation, molecular interactions which commonly appear both in ground and excited states, such as a part of dispersion interactions, do not appear in the formula.
ACS Paragon Plus Environment
10
Page 11 of 45 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
Figure 2. Excitations in the first-order interacting space and their interpretations.
To include interactions, we extended the excited-state wave function by including the excitation operators defined by the FOIS. We start from the CIS solution for the (QM+sQM)/MM system because CIS is formally an Ο ( N 4 ) method. Hence, the CIS computation including the sQM region is practically available. In addition, the accuracy of the variation solution for the single excitation manifold of the FOIS is better than that of the perturbation solution. The wave function to be solved becomes
ex Ψ iv-vi = CIS +
solvent
solute
∑ ∑ ∑ Sˆ
bj
A
A Sˆai HF dbjai
. (5)
b , j∈ A a , i
A This equation includes the interactions (iv)-(vi). The coefficients {dbjai } were solved by a second-
order perturbation theory using the CIS state for the zero-th order problem. To avoid the singularity problem in the double excitation manifold, we did not use an internally-contracted form of the wave function. The energy correction given by CIS and FOISPT2 calculations, ∆ CIS+PT2 = ∆ CIS + ∆ PT2
, (6)
ACS Paragon Plus Environment
11
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 12 of 45
was used for the second term of r.h.s. of eq. (3). In addition, the energy correction is easily decomposed into contribution from each solvent molecule because MOs are localized within each solvent and because the second order perturbation theory was adopted for solving the equations. To decompose the PT2 energy into the electronic polarization and dispersion effect in the excited state, we first solved the wave function, solvent Ψivex = 1 + ∑ ∑ Sˆbj dbjA CIS A b , j∈ A
(7)
to obtain the electronic polarization energy,
∆ iv
dispersion effect was estimated by subtracting
, with second-order perturbation theory. The ∆ iv
from ∆ PT2 . These terms could also be
decomposed into contributions from each solvent molecule. The limitation of the present definition is that the decomposition depends on the definition of fragments using LMOs. As both dispersion and polarization interactions originate from double excitations, calculated interaction energies for dispersion and polarization are dependent on the MOs and, therefore, interconvertible to each other. In addition, because our polarization energy was determined at the 2-nd order perturbation level, the electronic structure of the solute is the same as its ground state. Methods such as EFP21 and MMpol18 can handle polarization effect in the excited states in a self-consistent manner. To distinguish with complete polarization, the present definition would be precisely termed as first-order polarization effect.
2.2 Accelerated integral transformation code for the FOISPT2 calculation
ACS Paragon Plus Environment
12
Page 13 of 45 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
A computational program has also been extensively developed for applications to much larger systems. The most time-consuming part of the calculation was integral transformation from atomic-orbital basis to molecular-orbital basis, which limits the target of the application to small solute-solvent clusters (up to 20 solvent molecules) in our previous study.48 In the present study, an OpenMP-parallelized code for the atomic-orbital direct four-index transformation was implemented into the SMASH program,53 which is highly efficient and parallelized code for calculating two-electron repulsion integral over atomic orbitals. We also significantly reduced disk usage in the new code. These implementations realized fast calculations for realistic computational models. For example, one of the largest calculations in this study is a deprotonated retinal Schiff base in MeOH environment over 3700 basis sets (more than solvent molecules within 7 Å from the chromophore are included in the sQM part). This calculation completed in 5 days with Xeon (E5 2690 v3 2.6 GHz) 24 cores. With our old code, estimated computation time is around 244 days. After the CIS calculations with the LMO basis, a four-index transformation of the twoelectron repulsion integrals was performed to set up the molecular integrals necessary for the PT2 calculation. Because FOISPT2 calculates excited-state molecular interactions between fragments and chromophore, limited types of molecular integrals appear in the formula for the second-order perturbation equation: when a CIS state is used for the reference state of PT2 calculation, molecular integrals of (OIOI|O0V0) and (VIVI|V0O0) types become necessary. Here, “O” and “V” refer to occupied and vacant orbitals, respectively, while superscripts I and 0 denote fragment I and chromophore (fragment 0), respectively. For particular sets of atomic orbitals, the repulsion integrals were evaluated on-the-fly and the first indices are transformed:
( χ r χ s | χt χu ) → (φi χ s | χt χu ) , (φa χ s | χt χu ) , where {χ} and {φ} are AOs and MOs, respectively.
ACS Paragon Plus Environment
13
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 14 of 45
These transformed integrals are stored on disk for the next transformation, while AO integrals,
( χ r χ s | χt χu ) , are discarded to save disk space. The transformations of the second to fourth indices were accelerated by the DGEMM subroutine in the BLAS library.
2.3 MO localization MO localization scheme employed here was reported in our previous study,47, 54 and a brief description is given in this subsection. Above mentioned calculation is based on the fragmentation of total system into subsystems such as amino acid residues, peptide bond units, and solvent molecules as shown in Figure 1b. To localize MOs into these predefined fragments, localization schemes that use internal indices55-56 are not applicable. In our localization scheme, MOs were transformed to maximize overlap with externally-introduced reference MOs of each fragment using the minimum orbital deformation method.57. AO MO
M ij = ∑∑ φiA χ s CsCMO (8) , k Wk , j s
k
Here, φiA , a reference MO localized within a fragment A, is obtained by HF calculation of fragment A as an isolated system. { χ s } and {CsCMO , k } are AO and canonical MO coefficients of total system., respectively. {Wk , j } is a transformation matrix to be determined. For an amino acid residue fragments and a solvent molecule fragment, Hartree-Fock canonical MOs were calculated as an isolated species and adopted for reference molecular orbitals. In Figure 1b, a simple model peptide is shown for explanation. For fragments 1-3, reference MOs are prepared. A transformation matrix was determined by maximizing overlap integrals between LMOs and reference MOs. For the other fragment, such as peptide bond (-CO-NH-) units, MOs were
ACS Paragon Plus Environment
14
Page 15 of 45 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
localized using the Pipek-Mezey (PM) method58. Because the PM method can separate π orbitals from σ orbitals, the LMO for the peptide bond unit looks very similar to HF canonical MOs.
3. Computational details 3.1 para-Coumaric Acid in neutral and anionic forms The structures of p-CA in s-cis-anti and s-cis-syn conformers, [p-CAc]- (carboxylate form), and [p-CAp]- (phenolate form) are shown in Figure 3. A classical molecular dynamics (MD) calculation was performed to sample the structure of p-CA and [p-CAc]- in water solution. For the initial structure, p-CA was optimized in the gas phase. Periodic boundary conditions with a rectangular cell of 49.750 × 44.146 × 40.862 Å3 and 54.049 × 50.018 × 46.150 Å3 were set for pCA and [p-CAc]-, respectively. Moreover, a total of 2214 and 4273 TIP3P water molecules were positioned around the solute in each simulation cell, respectively. For the [p-CAc]- system, Na+ was additionally included to neutralize the system. The particle-mesh Ewald method was used to compute long-range electrostatic interactions. The general Amber force field59 were employed for p-CA and [p-CAc]-, while an Amber ff14SB force field60 was applied for water molecules. MD calculations were performed in an NPT ensemble at a constant pressure of 1 bar. The time step was 0.5 fs with a total simulation time of 5 ns. To equilibrate the system, the initial temperature of 500 K was reduced by an annealing method and maintained at the equilibrium temperature of 300 K monitored by a Nosé–Hoover thermostat. A total of 4 configurations were randomly acquired from the obtained trajectory (3–5 ns). Although statistics from just four snapshots is not meaningful, averaged values were calculated for comparing with the experimental data and for the discussion of the trend in the result. The results of each snapshot
ACS Paragon Plus Environment
15
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 16 of 45
are also given in Table below. For the calculations of the QM and sQM regions, the 6-31G* basis sets were used for p-CA and [p-CAc]-, and the 6-31G sets were employed for water molecules. Amber1461 and Gaussian0962 program packages were utilized for the classical MD simulation and geometry optimization, respectively.
Figure 3. Structure of p-CA in (1) s-cis-anti and (2) s-cis-syn forms, and (3) [p-CAp]-, and (4) [pCAc]-.
For each sampled configuration, QM/sQM/MM calculations were performed to compute the excitation energy. The QM region includes only chromophore, namely p-CA and [p-CAc]-. To define the sQM region, we used a parameter rsQM , the distance from the chromophore. The water molecule was included into the sQM region if at least one of its atom was within rsQM Å from the chromophore. The parameter rsQM was extended up to 7 Å. In our previous study, most of the first-order molecular interactions were recovered for an sQM region defined with rsQM = 7 Å.46 For the MM region, we considered water molecules within 30 Å from the chromophore. For the [p-CAc]- case, Na+ was found in a region of rsQM = 7 - 30 Å and included in the MM region.
ACS Paragon Plus Environment
16
Page 17 of 45 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
3.2 Retinal Schiff base The atomic coordinates of deprotonated retinal Schiff base (DPSB) in MeOH solution and protonated retinal Schiff base (PSB) in bR protein were obtained from our previous studies.46, 63 For the bR case, the initial structure was taken from the X-ray structure (PDB ID: 1C3W).64 The QM segment includes a retinal Schiff base with the side chains of Lys216, Asp85 residue modeled by CH3COO- (counterion group), and a proximal water molecule with a hydrogenbonding to Asp85; see active site (AS) model in Figure 4a. The QM/MM boundary was set on the Cβ-Cγ bond of Lys296 and on the Cα-Cβ bond of Asp85. Additionally, the hydrogen atom was used for a link atom. The atomic coordinates for the whole bR protein were optimized with QM(B3LYP)/MM(AMBER99) calculation. For the quantum chemical calculations, the 6-31G* and 6-31G basis sets were applied for the atoms in the AS region (see Figure 4a) and other atoms, respectively. For the MeOH solution, the four initial structures were taken from MD trajectory; for details, see reference.46 Although statistics from just four snapshots is not meaningful, averaged values were calculated for comparing with the experimental data and for the discussion of the trend in the result. The results with each snapshot are also given in Table below. The MeOH molecules located up to 20 Å away from the DPSB were considered in the model. The QM and MM regions consist of DPSB and MeOH molecules, respectively. To relax the structure of DPSB, QM(B3LYP)/MM(AMBER99) optimization was performed. The positions of the MM atoms were frozen during the optimization to retain its solvation structures. For the quantum chemical calculations, the 6-31G* and 6-31G basis sets were used for DPSB and MeOH, respectively.
ACS Paragon Plus Environment
17
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 18 of 45
(a) AS model
(b) MeOH-3 Å model
(c)
MeOH-7
Å
(d) bR-3 Å model
(e) bR-7 Å model
Figure 4. Computational models for a retinal Schiff base in MeOH solution and in bR protein environment. (a) Active site (AS) model (PSB with a water molecule and an acetate anion), (b) MeOH-3 Å model, (c) MeOH-7 Å model, (d) bR-3 Å model, and (e) bR-7 Å model. See main text for details.
To investigate the quantum mechanical effect on the calculated excitation energy, the sQM region was defined. MeOH-3Å and MeOH-7Å models were constructed for the MeOH environment. These models include MeOH molecules within 3 and 7 Å from DPSB in the sQM region; see Figures 4b and 4c. A total of 32 MeOH molecules were included in the sQM region for the MeOH-3Å model, while 126-136 MeOH molecules were included for the MeOH-7Å model (the number of MeOH molecules depends on the configurations employed). The rest of
ACS Paragon Plus Environment
18
Page 19 of 45 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 atoms in the model were included in the MM region. For the bR case, amino acid residues, peptide bonds, and water molecules that are within 3 and 7 Å from PSB were included in the sQM region for the bR-3Å and bR-7Å models, respectively, see Figures 4d and e. A total of 28 amino acid residues and 2 water molecules were included in the sQM region for the bR-3 Å model, while 64 residues plus 6 water molecules were in the bR-7Å models. To terminate the QM region, C and N terminals in the peptide bond were replaced by a hydrogen atom with a CαH bond length of 1.09 Å. The basis sets 6-31G(d)65-66 and 6-31G
were used for DPBS and PBS, and methanol
molecules and amino acids, respectively. An in-house modified version of the Gaussian 0962 was used for the CIS calculations with non-canonical HF MOs.
4. Results and Discussion 4.1 para-coumaric acid in water solution A solvatochromic molecule p-CA is a photoreceptor in photoactive yellow protein. Depending on the solvents and pH, the absorption maximum shifts about 150 nm.67-68 Therefore, it is one of the ideal molecules to investigate for the applicability of our computational method to solvatochromic systems. For the neutral species, the excited states of p-CA in gas phase and aqueous solution were investigated. For the anion species, [p-CAp]- in gas phase and [p-CAc]-, in aqueous solution were chosen because [p-CAp]- and [p-CAc]- dominantly appears in gas phase and aqueous solution. For p-CA in gas-phase, SAC-CI calculations were performed using the B3LYP optimized geometry. According to a previous study,69 the ground-state energies of four conformers of p-CA
ACS Paragon Plus Environment
19
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 20 of 45
are very close to each other (within 1 kcal/mol). The most stable one was the s-cis-anti conformer. The s-cis-syn conformer comes second with only 0.09 kcal/mol higher energy than the s-cis-anti conformer; see Figure 3 for the structure of the two conformers. In our B3LYP/631G* calculations, similar results were obtained; the s-cis-anti conformer is slightly more stable than the s-cis-syn (∆G = 0.6 kcal/mol at 298 K). For the π–π* excited state with the greatest oscillator strength, the excitation energies calculated with the SAC-CI method were 4.42 eV and 4.40 eV for the s-cis-syn and s-cis-anti conformers, respectively. These states were characterized as HOMO-LUMO excited state. In a previous CASPT2 study, excitation energy for the HOMOLUMO state was 4.49 eV and 4.52 eV for the s-cis-syn and s-cis-anti conformers, respectively.69 Experimentally observed absorption peak is at 4.04 eV70, and the SAC-CI result is overestimated by 0.3 eV. According to the CASPT2 study, this discrepancy could be improved by the basis sets.69 To evaluate the FOIS correction for the excited-state molecular interactions, CIS and PT2 calculations were performed. The CIS excitation energy, correction by CIS ( ∆ CIS ), and correction by CIS and PT2 ( ∆ CIS+PT2 ) are shown in Table 1. The CIS calculation with rsQM = 0 Å means that all the solvent molecules were described by the MM method. With rsQM = 7 Å, the CIS excitation energy was 5.02-5.11 eV (5.09 eV in average). As the ∆ CIS values show, the calculated excitation energy decreased by 0.07-0.23 eV (0.16 eV in average) relative to that with rsQM = 0 Å. With the CIS calculation, orbital delocalization and excitonic coupling effects were introduced into the calculated excitation energy. By considering the PT2 contribution, ∆ PT2 , the calculated excitation energy further decreased by 0.04-0.05 eV (0.04 eV in average). The PT2 correction originated from electronic polarization and dispersion effects. The average total shift, ∆ CIS+PT2 became -0.12
ACS Paragon Plus Environment
20
Page 21 of 45 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
eV to -0.28 eV (-0.20 eV in average). The result of the four snapshots indicates that the distribution of ∆ PT2 was less dependent than that of ∆ CIS . Figure 5 shows the spatial distribution of the PT2 contribution for snapshot 1 in Table 1. The dispersion effect dominantly contributed to the PT2 correction, while the polarization effect was negligible. This trend is common to all of the snapshots as shown in Figure S1a in Supporting information (SI). In addition, the largest contributions to the excitation energy arise from the solvent molecules in direct contact with the π electron system. After these FOISPT2 corrections were combined, the SAC-CI/AMBER + ∆ CIS+PT2 result was 3.96-4.25 eV (4.08 eV in average), which is in reasonable agreement with the experimental observation (4.00 eV)71 and the CASPT2 +ASEP/MD result (4.00 eV).69
ACS Paragon Plus Environment
21
Journal of Chemical Theory and Computation
Table 1. Result of the SAC-CI and FOISPT2 calculation of π-π* excited state of p-CA in gas phase and aqueous solution. Units are in eV.
gas snap shot
aqueous phase CIS/AMBER
SAC-CI
3Å
rsQM
rsQM
SAC-CI/ AMBER rsQM
7Å
7Å
0Å
∆
rsQM 0Å
7Å
∆
+∆
+∆
rsQM
rsQM
7Å
7Å
4.42g 4.40h
experiment / theory 4.06i / 4.49g, 4.19j, 3.94k
#1 #2 #3 #4
5.34 5.22 5.24 5.20
5.30 5.19 5.12 5.06
5.11 5.15 5.09 5.02
-0.23 -0.07 -0.15 -0.18
-0.28 -0.12 -0.19 -0.22
4.36 4.37 4.21 4.18
4.13 4.30 4.06 4.00
4.08 4.25 4.02 3.96
ave.e
5.25
5.17
5.09
-0.16
-0.20
4.28
4.12
4.08
4.00a / 4.00f ,4.03k
a
Experiment. Ref.[68]. bExcitation energy calculated with 7 Å minus that with 0 Å. CIS-AMBER level. c FOISPT2 correction was added to ∆ . d“∆ ”and “∆ ” denote that the calculated excitation energy with SAC-CI/AMBER ( 0 Å) was corrected by ∆ and ∆ , respectively. eAverage of four snapshots. fASEP/MD result with the polarization effect. CASPT2 for the electronic structure calculation. Ref.[69]. gs-cis-anti form. hscis-syn form. iExperiment. Ref.[70], jB3LYP with RISM-SCF-SEDD. Ref.[72]. kGreen’s Function theory. Ref. [73].
Indices of solvent molecules
Contribution to excitation energy / eV
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 22 of 45
4
22 22 Polarization
2
Dispersion
4
2
Figure 5 Decomposition analysis of the second-order interaction energy into each solvent’s contribution. Inset shows the position of water molecules, 2, 4, and 22, which contributed the greatest among the solvents. Result from snapshot 1 in Table 1.
ACS Paragon Plus Environment
22
Page 23 of 45 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 result of [p-CAp]- in gas phase is summarized in Table 2. Excitation energy by SAC-CI was 2.93 eV, which is comparable with that obtained with CASPT2 (2.89 eV)69 and experiment (2.88 eV).71
In aqueous solution, calculations were performed for [p-CAc]-, because the
carboxylate form is much more stable in aqueous phase. The SAC-CI/AMBER result distributes from 4.50 eV to 5.14 eV (4.81 eV in average), while the result was improved to 4.33-4.90 eV (4.56 eV in average) by including the FOISPT2 corrections. In experiment. absorption peak was observed at 4.35 eV.71 Regarding the corrections due to the interactions, ∆ CIS at rsQM = 7 Å distributes from -0.18 eV to -0.22 eV, while ∆ PT2 does from -0.04 eV to -0.05 eV. Similar to the neutral case, the distributions of ∆ CIS was much larger, and the ∆ PT2 correction was less dependent on structures. In Figure S1b of SI, the spatial distribution of the FOISPT2 corrections were presented. Similar to the p-CA case, the largest contributions to the excitation energy arise from the solvent molecules in contact with the π electron system.
ACS Paragon Plus Environment
23
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 24 of 45
Table 2. Result of the SAC-CI and FOISPT2 calculations of π-π* excited state of [p-CAp]- in gas phase and [p-CAc]- aqueous solution. Units are in eV. gas
snap shot
CIS/AMBER SAC-CI
rsQM 0Å
7Å
aqueous solution ∆ ∆
SAC-CI/ AMBER
+∆
+∆
rsQM
rsQM
rsQM
rsQM
rsQM
7Å
7Å
0Å
7Å
7Å
experiment/ theory 2.88a / 2.89h
2.93 #1 #2 #3 #4
5.41 5.60 5.38 5.43
5.19 5.42 5.16 5.23
-0.22 -0.18 -0.22 -0.20
-0.27 -0.24 -0.27 -0.25
4.60 5.14 4.62 4.89
4.38 4.96 4.40 4.69
4.33 4.90 4.35 4.64
ave.g
5.46
5.25
-0.21
-0.26
4.81
4.61
4.56
4.35c / 4.76h, 4.27i
a
Experiment. Reference [71]. cExperiment. Reference [67]. dExcitation energy calculated with 7 Å minus that with 0 Å. CIS-AMBER level. e FOISPT2 correction was added to ∆ . f “∆ ”and “∆ ” denote that the calculated excitation energy with SACCI/AMBER ( 0 Å) was corrected by ∆ and ∆ , respectively. gAverage of four snapshots. hASEP/MD result with the polarization effect. CASPT2 for the electronic structure calculation. Ref.[69]. iB3LYP with RISM-SCF-SEDD. Ref.[72].
ACS Paragon Plus Environment
24
Page 25 of 45 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
4.2. Solvatochromic shift of retinal chromophore Opsin shift is a solvatochromic shift of retinal Schiff base from MeOH solution to protein environments.74 The experimentally observed absorption energy in MeOH is, 3.41 eV (237 nm),75 while those in protein environments spread over a wide energy range. One of the most significant bathochromic shifts is observed in bR protein environment (2.18 eV, 569 nm),76 with a shift of 1.23 eV. In this study, we applied QM/sQM/MM method to this solvatochromic shift. 4.2.1 The rsQM dependence of the sQM effect Table 3 shows the excitation energies of the first singlet excited state of DPSB and PSB, respectively. Their rsQM dependence is also summarized. CIS and its PT2 correction were used to account for the sQM effect. As rsQM increases, calculated excitation energy decreases. At rsQM = 7 Å, the change in the excitation energy became very small. This trend is commonly observed in our previous calculations.46-47 Compared with the result when rsQM = 0 Å (no sQM region), calculated excitation energy decreased by 0.07-0.12 eV (0.10 eV in average) in MeOH case and by 0.17 eV in bR case. The PT2 correction also reduced the excitation energy by 0.07-0.10 eV (0.09 eV in average) in MeOH case and by 0.18 eV in bR case. The total CIS + PT2 corrections, ∆ CIS+PT2 ,
were reduced by 0.14 -0.21 eV (0.18 eV in average) in MeOH case and by 0.35 eV in bR
case. These contributions are significant, not only in the excitation energy, but also in the relative change (-0.17 eV) from the MeOH solution to the bR protein environment. These corrections are added to the SAC-CI correlated calculation using equation 3 as the sQM contribution to the excitation energy.
ACS Paragon Plus Environment
25
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 26 of 45
Table 3. Result of the SAC-CI and FOISPT2 calculation of π-π* excited state of DPSB and PSB in gas phase, MeOH solution, and bR protein environment. Units are in eV. Opsin shift at several sQM values are also shown. Gas snap shot
in solution / protein environment CIS/AMBER
∆
∆
SAC-CI/ AMBER
∆
rsQM
rsQM
rsQM
rsQM
rsQM
rsQM
SAC-CI
0Å (1) DPSB in MeOH solution 3.52 #1 4.18 #2 4.17 #3 4.19 #4 4.24 ave.e 4.20 (2) PSB in bR protein 3.39 (3) Opsin shift
Exptl.
3Å
7Å
3Å
7Å
3Å
7Å
0Å
7Å
7Å
4.13 4.08 4.11 4.15 4.12
4.11 4.06 4.10 4.12 4.10
-0.05 -0.09 -0.08 -0.09 -0.08
-0.07 -0.11 -0.09 -0.12 -0.10
-0.04 -0.08 -0.06 -0.07 -0.06
-0.07 -0.10 -0.08 -0.09 -0.09
3.49 3.47 3.46 3.54 3.47
3.42 3.36 3.37 3.42 3.37
3.36 3.26 3.29 3.33 3.28
3.41a
3.29
3.22
-0.10
-0.17
-0.17
-0.18
2.36
2.19
2.01
2.18f
1.27
1.23
1.11 a
∆
75
b
Experiment. Reference [ ]. Change in calculated excitation energy from 0 Å at CISAMBER level. c FOISPT2 correction. d “∆ ”and “∆ ” denote the SAC-CI/AMBER excitation energy corrected by ∆ and ∆ , respectively. eAverage of four snapshots. f Experiment. Reference [76].
ACS Paragon Plus Environment
26
Page 27 of 45 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
As a general trend, the environmental effect calculated with CIS requires a relatively large rsQM value because the excitonic coupling effect is introduced by the single excitations. As -3 indicated by the Förster theory77, the excitonic coupling interaction qualitatively scales as rsQM .
In the MeOH environment, the calculated CIS excitation energies with rsQM = 0 Å, 3 Å, and 7 Å were 4.20 eV, 4.12 eV, and 4.10 eV, respectively; average of snapshots, see Table 3. The major part of the change was introduced at rsQM = 3 Å. The amount of the change became small when rsQM was extended to 7 Å. In our previous studies, the change almost converged at rsQM = 6 Å.4647
However, in the bR environment, a larger rsQM value was necessary. By extending the rsQM
value from 3 Å to 5 Å, the calculated excitation energy further decreased by 0.09 eV. The result can be interpreted by the Förster theory.77 Because the excitonic coupling depends on the magnitude of the transition dipole moment, amino acid residues in the protein environment contributed to the excitation energy in such a long-range behavior. As shown in Table 4, the calculated excitation energy became stable at rsQM = 7 Å. In contrast, the PT2 correction showed a rapid saturation behavior with increasing rsQM . The result at rsQM = 3 Å was within the deviation of 0.01 eV (average of the snapshots) from that at rsQM = 7 Å. In the MeOH solution, the result of the MeOH-3Å model (-0.063 eV) amounts to 73
% of that from the MeOH-7Å model (-0.086 eV); average of four snapshot, see Table 3. A similar trend was observed in the protein environment as shown in Table 3. The PT2 energy in the bR-3 Å model was -0.17 eV, which already accounts for 94 % of that in the bR-7 Å model. As these behavior shows, the FOISPT2 corrections due to the environment are covered by the sQM region with rsQM = 7 Å for the systems in the present study.
ACS Paragon Plus Environment
27
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 28 of 45
Calculated and experimental excitation energies for the first excited state of DPSB in MeOH and PSB in bR protein were summarized in Table 3. For DPSB in MeOH solution, the calculated excitation energy with rsQM = 7 Å was 3.28 eV (average of 4 snapshots), while experimentally observed one is 3.41 eV.75 For PSB in bR protein environment, the calculated excitation energy with rsQM = 7 Å was 2.01 eV, while the experimentally observed one is 2.18 eV.76 In both models, the calculated excitation energies underestimated the experimental data by 0.13 eV and 0.17 eV for MeOH-7Å and bR-7Å models, respectively. However, compared with the QM/MM calculation without the sQM region, the calculated solvatochromic shift, relative excitation energy, showed an improvement. The calculated opsin shift at rsQM = 7 Å was 1.27 eV, while the experimental one is 1.23 eV.
4.2.2. Decomposition analysis of the second-order energy The PT2 energy was decomposed into polarization and dispersion contributions. As summarized in Figure 6, dispersion effect is much greater than polarization effect and decreased excitation energies in the MeOH-7Å and bR-7Å models by 0.09 eV and 0.17 eV, respectively. In contrast, the polarization effect decreased excitation energy only by 0.001 eV and 0.01 eV for the MeOH-7Å and bR-7Å models, respectively. Regarding the polarization effect, there are currently two different results from the previous computational studies. The first used a polarizable force field model for the protein environment and produced red shift of 0.34 eV30 and 0.14 eV.
45
On
the other hand, a polarized LoProp model for rhodopsin and a three-layer QM/QM/MM calculation gave red shifts of 0.06 and 0.03 eV for the electronic polarization.41 The result of our QM/sQM/MM calculation is closer to the latter, suggesting that electronic polarization effect has minor influence in the solvatochromism of this system.
ACS Paragon Plus Environment
28
Page 29 of 45 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
Figure 6. CT+Excitonic, dispersion, and polarization contributions to the excitation energies of DPSB in MeOH and PSB in bR. Units are in eV. Averaged values were used for DPSB in MeOH. Other theoretical estimate for polarization effect of PSB in bR was -0.34 eV30 and -0.03 eV.41
The dispersion and polarization energies were further decomposed into the contribution from each solvent molecule. The result for the MeOH-3Å model was investigated. As seen in Figure 7a, several solvent molecules presented large contributions. The shortest atom-to-atom distances between DPSB and each MeOH molecule were given in Figure 7b. In total, seven MeOH molecules are within 3 Å from DPSB, illustrating a clear correlation with the magnitude of the dispersion contribution. The amount of stabilization from these MeOH molecules was -0.057 eV, accounting for 68 % of the total dispersion energy in the MeOH-3Å model. Figure 7c shows the locations of the seven MeOH molecules. These MeOH molecules interact with the DPSB chromophore at the π-skeleton from the top and bottom, suggesting that this direct interaction caused a large dispersion energy.
ACS Paragon Plus Environment
29
Journal of Chemical Theory and Computation
(a) Energy decomposition to each solvent molecule Fragment index
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 0.000
Second-order contributions/eV
-0.002 -0.004 -0.006 -0.008 -0.010
Tot -0.084 eV Pol -0.0002 eV Disp -0.084 eV
-0.012 -0.014
(b) Distance from DPBS 7
Distance from π-skelton/Å
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 30 of 45
6 5 4 3 2 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
Fragment index
(c) MeOH molecules with large contributions (-0.004 eV)
Figure 7. Result of the decomposition analysis for the MeOH-3Å model. (a) Contributions of each MeOH molecule. (b) Distance from each MeOH molecule to DPSB. For the distance, the shortest atom-to-atom distance was used. (c) MeOH molecules with large contributions (< -0.004 eV).
ACS Paragon Plus Environment
30
Page 31 of 45 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 decomposition analysis was also performed for the bR case. Dispersion and polarization contributions were decomposed into each fragment (amino acids, peptide-bond units, and water molecules). The result is shown in Figure 8a. The dispersion contribution arose from several fragments, the same trend as the MeOH case. Five amino acids contributed large stabilization (< -0.01 eV). These gave a total stabilization of -0.12 eV, which accounts for 72 % of the total PT2 energy in the bR-3Å model. The largest and second largest contributions were given by 19-th (Tyr185) and 8-th (Trp86) fragments with contributions of -0.054 and -0.028 eV, respectively. The locations of the five amino acids were visualized in Figure 8b. Tyr185 exists just above the retinal chromophore with their molecular planes parallel to each other, indicating that their πelectrons are very closely located.
ACS Paragon Plus Environment
31
Journal of Chemical Theory and Computation
(a) Energy decomposition analysis to each fragment Fragment index 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 0.000 -0.010
Second-order contributions/eV
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 32 of 45
Trp182
Met118 Thr90
-0.020 -0.030
Trp86
-0.040 -0.050
Tot -0.17 eV Pol -0.0089 eV Disp -0.16 eV
Tyr185
-0.060
(b) Amino acids with large contributions (
