Novel Approach to Excited-State Calculations of Large Molecules

Article pubs.acs.org/JPCB

Novel Approach to Excited-State Calculations of Large Molecules Based on Divide-and-Conquer Method: Application to Photoactive Yellow Protein Takeshi Yoshikawa,† Masato Kobayashi,‡ Atsuhiko Fujii,† and Hiromi Nakai*,†,§,∥,⊥ †

Department of Chemistry and Biochemistry, School of Advanced Science and Engineering and §Research Institute for Science and Engineering, Waseda University, Tokyo 169-8555, Japan ‡ Waseda Institute for Advanced Study, Waseda University, Tokyo 169-8050, Japan ∥ CREST, Japan Science and Technology Agency, 4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan ⊥ ESICB, Kyoto University, Kyotodaigaku-Katsura, Nishigyoku, Kyoto 615-8520, Japan S Supporting Information *

ABSTRACT: In this study, the divide-and-conquer (DC) method is extended to configuration-interaction singles, time-dependent density functional, and symmetry-adapted cluster configuration interaction (SACCI) theories for enabling excited-state calculations of large systems. In DC-based excited-state theories, one subsystem is selected as the excitation subsystem and analyzed via excited-state calculations. Test calculations for formaldehyde in water and a conjugated aldehyde demonstrate the high accuracy and effectiveness of these methods. To demonstrate the efficiency of the method, we calculated the π−π* excited state of photoactive yellow protein (PYP). The numerical applications to PYP confirm that the DC-SACCI method significantly accelerates the excited-state calculations while maintaining high accuracy.

density-fragment interaction approach10 enables large-scale TDDFT calculations by describing the interactions between many quantum-mechanical fragments. Atomic-orbital (AO) based TDDFT11,12 enables linear-scaling computation via prescreening techniques and sparse matrix algebra. There are several approaches that utilize localized molecular orbitals (MOs), including fragment-localized MO TDDFT,13 localexcitation-approximation CIS and TDDFT,14 and the local CIS approach.15 Also, linear-scaling TDDFT with a local density matrix16,17 has been proposed. A giant SACCI scheme18 was proposed for the excited-state calculation of molecular crystals. We developed fragmentation-based linear-scaling methods for ground-state theories based on the divide-and-conquer (DC) method. The DC method for density functional theory was first proposed by Yang and co-workers.19,20 Unlike other fragmentation-based methods, the DC method requires no artificial prediction related to the positions of the spin and/or charge because the distribution of electrons in the system under consideration is uniformly settled by the common Fermi level; this enables accurate treatment of large π-conjugated systems. We assessed the DC method using Hartree−Fock (HF) calculations of closed-shell21−24 and open-shell25,26 systems,

1. INTRODUCTION Many biochemical systems, such as rhodopsin in vertebrate eyes and photoactive yellow protein (PYP) in purple bacteria, are photoactive. The absorption wavelength of a photoactive protein is affected by the photon-absorbing pigment as well as interactions between the pigment and its neighboring residues, which cause red- or blue-shifts from the excitation energy of the free pigment. Although computational predictions of the shifts in biomolecular systems are in demand, quantum chemical excited-state calculations of large biomolecules are difficult because of the high computational cost. Recently, excited-state calculations based on symmetryadapted cluster configuration interaction (SACCI) theory,1,2 which is one of the most accurate excited-state theories, successfully explained and predicted the mechanisms of photoactive biomolecular systems using an active-site model selected from the system.3,4 Although active-site models work reasonably well for many biomolecular systems, their calculations overlook interactions between the pigment and outside residues, which may be important for accurate evaluation of the excitation energy. Several fragmentation-based approaches that enable reduced computation for excited-state theories have been reported. For example, the fragment molecular orbital (FMO) method was applied to the configuration interaction singles (CIS)5−7 and time-dependent density functional theories (TDDFT).8,9 The © 2013 American Chemical Society

Received: February 21, 2013 Revised: April 8, 2013 Published: April 9, 2013 5565

dx.doi.org/10.1021/jp401819d | J. Phys. Chem. B 2013, 117, 5565−5573

The Journal of Physical Chemistry B

Article

In eq 3, Cαp and εαp represent the subsystem MO coefficient and orbital energy, respectively; these values are determined by solving the following Roothaan-Hall equation for subsystem α:

and applied it to static and dynamic (hyper)polarizability calculations.27−29 Furthermore, we extended the DC method to electron-correlation theories, specifically the second-order Møller−Plesset perturbation,30−34 coupled-cluster,35,36 and symmetry-adapted cluster (SAC) theories,37 with the assistance of energy density analysis (EDA).38 Energy gradient schemes for the DC method were also recently established.39,40 For biomolecules, we successfully elucidated interactions between viral peptides and their inhibitors.41,42 Although the application of the DC method to excited-state theories is desired for the accurate prediction of the excitation energies of biomolecules, it has not yet been attempted. In the present study, we extended the DC methodology to excited-state theories, specifically CIS, TDDFT, and SACCI. In these DC-based excited-state theories, the subsystem that contains the chromophore or pigment is first selected as the excitation subsystem (ES). Excited-state calculation of the ES yields the excitation energies, which include the effect of factors outside the ES through the electrostatic potential. The organization of this article is as follows: Section 2 presents the theoretical aspects of the DC-CIS, DC-TDDFT, and DCSAC/SACCI methods after a brief summary of the DC-SCF method; section 3 shows the numerical applications of the present scheme via calculations of formaldehyde in water, conjugated aldehyde, and PYP; and section 4 comprises the conclusion.

Fα Cαp = εpα Sα Cαp

Here, Sα and Fα represent the local overlap and Fock matrices for subsystem α, i.e., submatrices of the entire overlap and Fock matrices in the basis of L(α). The Fermi level can be determined via constraint of the total number of electrons, ne, as follows: ne = Tr[DDCS] =

Dμν ≈

=

∑

(1)

n) ES diES( , a |Φi , a ⟩

a

(7)

ES ES ̂ AiaES, jb = ⟨ΦiES , a |H − E0 |Φ j , b⟩

= (εaES − εiES)δijδab − ⟨j ES a ES || i ESbES⟩

(8)

ES(n) ES ̂ ES where EES , can 0 = ⟨Φ0 |H|Φ0 ⟩. The excitation energies, ω be obtained from the eigenvalues of this matrix. The dual-buffer DC scheme, which was first proposed in the post-HF electron-correlation calculation33 can be applied to the excited-state calculation. In this scheme, the Fock matrix, F, is first constructed from the density matrix converged with the preceding conventional HF calculation. Then, the RoothaanHall equation for the ES is solved to construct the MOs used for the excited-state calculation, as follows:

(2)

(3)

where Pα is the partition matrix with elements of

α Pμν

(6)

Here, |ΦES 0 ⟩ is defined as the Slater determinant of the occupied MOs in the ES that have lower energies than the Fermi level determined via the DC-HF calculation; the subscripts {i, j} and {a, b} refer to occupied and virtual MOs, respectively. Hereafter, dES(n) is referred to as the CI coefficients that are deduced as the normalized eigenvectors of the Hamiltonian matrix, AES, which is expressed using the twoelectron integral notation of ⟨ij|ab⟩ = ∫ ∫ ϕi(r1)ϕj(r2)r−1 12 ϕa(r1)ϕb(r2)dr1dr2, as follows:

α α Dμν = Pμν ∑ fβ (εF − εpα)CμαpCναp

⎧1 [μ ∈ S(α) ∧ ν ∈ S(α)] ⎪ ⎪ = ⎨1/2 [μ ∈ B(α) ∧ ν ∈ S(α)] ⎪ ∨ [μ ∈ S(α) ∧ ν ∈ B(α)] ⎪ ⎩0 otherwise

∑ ∑ i

In this equation, Dα represents the local density matrix for subsystem α, which is obtained using the Fermi level, εF, and Fermi function, fβ(x) = [1 + exp(−βx)]−1, with an inverse temperature parameter, β, as follows: p

μ ∈ L(α)

occ(ES) vir(ES) ES(n) |ΨCIS ⟩=

α Dμν

α

(Dα Sα )μμ

Then, the entire density matrix, DDC, can be obtained from eq 2. The same procedure is adopted for the DC-DFT calculation by substituting the Fock matrix with the Kohn− Sham Hamiltonian. 2.2. DC-Based Excited-State Theory. The DC-based excited-state theory utilizes the subsystem MOs obtained from the preceding DC-HF or DC-DFT calculations. Following the ground-state DC calculation, the subsystem that includes the pigment is selected as the ES. The excited-state wave functions of the ES are constructed from its MOs, which can be classified as either occupied or virtual via the Fermi level. 2.2.1. DC−CIS/TDDFT. In the DC−CIS method, the nth excited-state wave function in the ES is expressed as the linear combination of singly excited determinants, |ΦES i,a ⟩, from the reference determinant, |ΦES 0 ⟩, as follows:

In the DC-HF method, the density matrix of the entire system, D, is constructed from the subsystem contributions, as follows: DC Dμν

∑ ∑ α

2. THEORY 2.1. DC-HF/DFT. In the DC method, the system under consideration is spatially divided into disjoint subsystems, which are called the central regions. A set of AOs corresponding to a central region, α, is denoted by S(α). To improve the description of the subsystem, the region neighboring the central region, which is called the buffer region, is considered during the expansion of the subsystem MOs. A set of AOs corresponding to the buffer region of subsystem α, which is denoted by B(α), is added to S(α) to construct a set of AOs in the localization region of subsystem α, L(α), as follows: L(α) ≡ S(α) ∪ B(α)

(5)

ES ES ES F ESCES p = εp S C p

(9)

The Fermi level that separates the occupied and virtual MOs is determined as the middle of the highest occupied and lowest virtual MO energies obtained from the preceding conventional HF calculation. Alternatively, it can also be defined using eq 6 by calculating the MOs in all the subsystems.

(4) 5566



Article ES

This procedure can also be applied to TDDFT calculations by substituting the CIS Hamiltonian of eq 8 with that of TDDFT.43 2.2.2. DC-SAC/SACCI. In a previous study,37 we proposed the ground-state DC-SAC theory. Here, we briefly summarize the conventional SAC/SACCI and DC-SAC theories and present the DC-SACCI method. The SAC expansion for a totally symmetric singlet ground state is expressed as |ΨSAC⟩ = exp(S)̂ |Φ0⟩

ES |ΨSAC ⟩ = exp(S ̂ )|Φ0ES⟩

where S represents the symmetry-adapted single and double excitation operator from the reference determinant, |ΦES 0 ⟩, in the ES, as follows: ES

Ŝ

+

(10)

occ vir

a

i,j

∑ ∑

(11)

EES SAC

Sî , a = (aa†↑ai ↑ + aa†↓ai ↓)/ 2

+ ab†↓ai ↓aa†↑aj ↑)/4

ES ES ES ES = E0ES + ⟨i ESj ES |a ESbES⟩(2ciES , a c j , b − ci , b c j , a + cij , ab

(13)

3 c′ijES, ab )

+

(24)

+ aa†↓ai ↓ab†↑aj ↑ − ab†↑ai ↑aa†↓aj ↓ − ab†↓ai ↓aa†↑aj ↑))/4 3

Based on the SAC wave function in the ES, the DC-SACCI expansion for the nth excited state is

(14)

Here, apσ and a†pσ represent the annihilation and creation operators, respectively, for spatial orbital p with σ spin. The SAC expansion coefficients, i.e., ci,a, cij,ab, and cij,ab ′ , are determined by solving the following nonvariational equations ⟨Φ0|Ĥ − ESAC|ΨSAC⟩ = 0

(15)

⟨Φ0|SÎ (Ĥ − ESAC)|ΨSAC⟩ = 0

(16)

ES(n) |ΨSACCI ⟩=

∑ dKES(n)P ̂

ES

ES †

ES R̂K |ΨSAC ⟩

K

(25)

{R̂ ES† K }

where represents a set of excitations, each of which generates a basis for the target state in the ES. P̂ ES = 1 − ES |ΨES SAC⟩⟨ΨSAC| is the projection operator that projects out the ground-state wave function in the ES. The DC-SACCI coefficients, dES(n) , and excitation energies, ωES(n), are calculated K by solving the following nonvariational equation:

̂ . Here ESAC is the SAC energy, which with SÎ = Ŝi,a, Ŝij,ab, or S′ij,ab comprises the HF energy, EHF, and the SAC correlation energy, ΔESAC. The SAC correlation energy can be calculated as

ES

ES ES(n) ⟨Φ0ES|R̂K (Ĥ − ESAC − ωES(n))|ΨSACCI ⟩=0

3 c′ij , ab )

(26)

3. RESULTS AND DISCUSSION In this section, the performance of the present DC-CIS, DCTDDFT, and DC-SACCI methods are numerically assessed by comparing the results with those of conventional CIS, TDDFT, and SACCI calculations. All CIS and TDDFT calculations were performed using a modified version of the GAMESS program,44 while SACCI calculations were performed using a modified version of the Gaussian09 program.45 In this section, only the excited-state calculations were performed using the DC scheme following conventional calculation of the ground states. All chemical core orbitals are frozen during the excited-state calculations. For the TDDFT calculations, the long-range corrected Becke−Lee−Yang−Parr (LCBLYP) exchange-correlation functional46−48 was adopted unless otherwise noted. The perturbation selection with a LevelThree threshold in the Gaussian program was used for the SACCI calculations. 3.1. Formaldehyde in Water. We first assessed the DCbased excited-state theory through calculations of the n−π* and

(17)

Based on the SAC wave function, the SACCI expansion for the nth excited state is expressed as †

(18)

where {R̂ †K} represents a set of excitations, each of which generates a basis for the target state. P̂ = 1 − |ΨSAC⟩⟨ΨSAC| is the projection operator that projects out the ground-state wave function. The SACCI coefficients, d(n) K , and the excitation energies, ω(n), are calculated by solving the nonvariational equation, as follows: n) ⟨Φ0|R̂K (Ĥ − ESAC − ω(n))|Ψ(SACCI ⟩=0

(23)

EES SAC

ES ES ESAC =E0ES + ΔESAC

S′̂ ij , ab = (2aa†↑ai ↑ab†↑aj ↑ + 2aa†↓ai ↓ab†↓aj ↓ + aa†↑ai ↑ab†↓aj ↓

K

(22)

Here, is the SAC energy of the ES. Note that is not the ground-state DC-SAC energy that was presented in the previous paper;37 instead, it is defined as

(12)

Siĵ , ab = (aa†↑ai ↑ab†↓aj ↓ + aa†↓ai ↓ab†↑aj ↑ + ab†↑ai ↑aa†↓aj ↓

∑ dK(n)PR̂ ̂K |ΨSAC⟩

(21)

a,b

ES ES ES ⟨Φ0ES|SÎ (Ĥ − ESAC )|ΨSAC ⟩=0

where

n) |Ψ(SACCI ⟩=

ES (cijES, abSiĵ , ab + c′ijES, ab S′̂ ijES, ab )

ES ES ⟨Φ0ES|Ĥ − ESAC |ΨSAC ⟩=0

a,b

ΔESAC = ⟨ij|ab⟩(2ci , acj , b − ci , bcj , a + cij , ab +

ES

̂ ciES , a Si , a

̂ES ̂′ES are the same as in eqs 12−14, but are defined ŜES i,a , Sij,ab, and Sij,ab using the annihilation and creation operators of the MOs in the ES† ES, i.e., aES pσ and apσ . The expansion coefficients are determined by solving the nonvariational equation for the ES, as follows:

∑ ∑ ci ,aSî ,a + ∑ ∑ (cij , abSiĵ , ab + c′ij , ab S′̂ ij , ab ) i

∑ ∑

i,j

S ̂ = SŜ + SD̂ occ vir

occ(ES) vir(ES)

=

i a occ(ES) vir(ES)

where Ŝ represents the symmetry-adapted single and double excitation operators, as follows:

=

(20)

E ̂S

(19)

In the DC-SAC/SACCI theory, the ground-state SAC wave function for the ES is 5567



Article

σ−π* excitation energies of a formaldehyde molecule in 16 explicit water molecules, i.e., H2CO + 16H2O, for which the coordinates are provided in Table S1 (Supporting Information). The ES comprises the H2CO molecule and atoms in the union of spheres centered at the four atoms of the H2CO molecule within radius rex. Figure 1 and Table 1 show the ES-

the deviations of the calculated excitation energies from experimental values are presented in parentheses. Empirically, the n−π* excitation energy of formaldehyde in the gas phase is 4.07 eV,50 while that in water is 4.28 eV;51 this indicates a blueshift of 0.21 eV. In the gas-phase calculations, the SACCI method gives the best result with a deviation of 0.12 eV. The deviation of the results of the TDDFT method is also small (0.15 eV), while the deviation of the results of the CIS method is relatively large. By increasing the number of water molecules, the excitation energies determined using the CIS and SACCI methods decrease gradually and approach the experimental value, while that determined using TDDFT increases. The SACCI result for H2CO + 61H2O gives the closest excitation energy (4.41 eV) to the empirical value with a shift from the gas-phase results of 0.22 eV. 3.2. Conjugated Aldehyde. Next, we applied the present methods to the n−π* excited-state calculations of a conjugated aldehyde, i.e., C16H17CHO, which is depicted in Figure 2. The CC, C−C, CO, and C−H lengths were fixed at 1.36, 1.46, 1.22, and 1.10 Å, respectively, and all atoms lie in the same plane. The ES consists of the terminal CHO group and several C2H2 units, the number of which is denoted by nex (see Figure 2, which shows an example for nex = 2). Figure 3 shows the ESsize (nex) dependence of the n−π* excitation energies obtained via DC-CIS, TDDFT, and SACCI calculations of C16H17CHO using the 6-31G** basis set. The conventional CIS, TDDFT, and SACCI excitation energies are represented by dashed lines in Figure 3. For all three methods, the results of the DC-based excited-state calculations show smooth convergence with respect to the ES size and deviate by an absolute value of 0.12 eV or less for the ES size of nex ≥ 3. The efficiencies of the DC-CIS and TDDFT methods were determined by measuring the central processing unit (CPU) time required using an Intel Xeon X5470 (3.33 GHz) processor on a single core. For diagonalization of the Hamiltonian matrix, the Davidson algorithm was adopted; the maximum number of expansion vectors used by the solver’s iterations was set at 50. The lowest 10 singlet excited states were obtained. Table 3 shows the CPU times for the DC-CIS and TDDFT calculations of C16H17CHO using the 6-31G** basis set. The times required for the preceding HF or DFT calculation are not included. It was confirmed that the DC method drastically reduces the computational time required with only a slight energy error. Table 4 shows the basis-set dependence of the n−π* excitation energies obtained via the DC-CIS, TDDFT, and SACCI calculations. 6-31G,52 6-311G,53 6-311G**,53 ccpVDZ,54 and cc-pVTZ54 basis sets were adopted in addition

Figure 1. ES-size (rex) dependence of the excitation energies obtained via DC-CIS, TDDFT, and SACCI calculations of the H2CO + 16H2O system using the 6-31G** basis set. The conventional CIS, TDDFT, and SACCI excitation energies are indicated by dashed lines.

size (rex) dependence of the n−π* and σ−π* excitation energies obtained via DC- and conventional CIS, TDDFT, and SACCI calculations of H2CO + 16H2O using the 6-31G** basis set.49 The conventional CIS, TDDFT, and SACCI excitation energies are shown as dashed lines in Figure 1. In Table 1, deviations of the DC excitation energies from those derived from conventional methods are presented in parentheses. For all three methods, the DC-based excitedstate calculations showed 0.16 eV or less absolute-value excitation-energy deviations for the ES size of rex ≥ 3.0 Å. The excitation energy deviations converge to zero when the ES is enlarged in the DC method, while the n−π* excitation energy deviation for the DC-TDDFT calculation with rex = 2.0 Å is smaller than that with rex = 2.5 Å; this is probably due to error cancellation. Next, we increased the number of water molecules in the calculation to investigate convergence with respect to the size of the solvation sphere. Here, the ES size was fixed at rex = 3.5 Å, and the 6-31G** basis set was adopted. Table 2 shows the dependence of the n−π* excitation energy of H2CO on the number of explicit water molecules; the gas-phase results calculated using conventional methods are also included, and

Table 1. ES-Size (rex) Dependence of the n−π * and σ−π* Excitation Energies (in eV) Obtained via DC-CIS, TDDFT, and SACCI Calculations of a H2CO + 16H2O System Using the 6-31G** Basis Seta CIS

TDDFT

SACCI

rex [Å]

n−π*

(diff.)

σ−π*

(diff.)

n−π*

(diff.)

σ−π*

(diff.)

n−π*

(diff.)

σ−π*

(diff.)

2.0 2.5 3.0 3.5 4.0 Conv.

3.51 4.64 5.01 5.09 5.09 5.10

(−1.58) (−0.46) (−0.09) (−0.00) (−0.01)

8.34 8.55 9.92 9.92 9.82 9.78

(−1.44) (−1.23) (+0.15) (+0.15) (+0.04)

3.00 2.84 4.17 4.23 4.23 4.24

(−1.24) (−1.40) (−0.07) (−0.01) (−0.01)

7.93 7.72 9.38 9.40 9.41 9.42

(−1.49) (−1.70) (−0.05) (−0.03) (−0.01)

4.31 4.30 4.52 4.55 4.61 4.68

(−0.37) (−0.38) (−0.16) (−0.13) (−0.07)

9.65 9.65 9.58 9.53 9.54 9.50

(+0.15) (+0.15) (+0.08) (+0.03) (+0.04)

a

The results of conventional CIS, TDDFT, and SACCI methods are shown at the bottom. The deviations of the calculated excitation energies from the conventional results are shown in parentheses. 5568



Article

Table 2. Dependence of the n−π* Excitation Energy (in eV) of H2CO on the Number of Explicit Water Molecules (n)a H2CO + nH2O CIS TDDFT SACCI Exptl.

n=0

(diff.)

n = 16

(diff.)

n = 20

(diff.)

n = 39

(diff.)

n = 61

(diff.)

4.58 3.92 4.19 4.07b

(+0.51) (−0.15) (+0.12)

5.09 4.23 4.55 4.28c

(+0.81) (−0.05) (+0.27)

5.03 4.43 4.46 4.28c

(+0.75) (+0.15) (+0.18)

5.01 4.46 4.45 4.28c

(+0.73) (+0.18) (+0.17)

4.92 4.61 4.41 4.28c

(+0.65) (+0.33) (+0.13)

a

The gas-phase results (n = 0) obtained via conventional calculations are also presented. The deviations of the calculated excitation energies from experimental values are presented in parentheses. bref 50. cref 51.

to the 6-31G** set; note that conventional SACCI with ccpVTZ could not be performed due to the large computational cost. The ES size was fixed at nex = 3. In Table 4, the deviations of the DC excitation energies from those obtained using conventional methods are presented in parentheses. The energy deviations do not show a significant dependence on the adopted basis set and, at a maximum of 0.15 eV, are quite small. Therefore, the use of a larger basis set does not reduce the effectiveness of the present method. Finally, the accuracies of the DC-TDDFT excitation energies obtained using typical functionals, i.e., SVWN,55,56 BLYP,46,47 B3LYP,57 and LCBLYP, were examined. Table 5 shows the

Figure 2. Structure of the C16H17CHO molecule and schematic of the ES in the DC calculation.

Table 5. DFT Functional Dependence of the n−π* Excitation Energy (in eV) of a Conjugated Aldehyde C16H17CHO Obtained Using DC- and Conventional TDDFT Methodsa

Figure 3. ES-size (nex) dependence of the excitation energies obtained via DC-CIS, TDDFT, and SACCI calculations of C16H17CHO using the 6-31G** basis set. The conventional CIS, TDDFT, and SACCI excitation energies are indicated by dashed lines.

CIS

TDDFT

1 2 3 4 Conv.

14 49 179 242 708

114 402 1290 1944 4530

conv.

DC

(diff.)

SVWN BLYP B3LYP LCBLYP HF

2.15 2.31 3.17 3.53 4.48

2.30 2.45 3.15 3.50 4.42

(+0.15) (+0.14) (−0.03) (−0.03) (−0.06)

a

The deviations of the energies from the conventional results are shown in parentheses. An ES size of nex = 3 was adopted for the DC calculations.

Table 3. CPU Times (in s) of the DC-CIS and TDDFT Calculations of C16H17CHO Using the 6-31G** Basis Seta nex

functional

functional dependence of the n−π* excitation energies obtained via DC (with nex = 3) and conventional TDDFT methods. The excitation energies obtained using pure DFT functionals (i.e., SVWN and BLYP) were lower than those obtained using wave function-based correlated SACCI, i.e., 3.87 and 3.76 eV by the conventional and DC-SACCI methods, respectively (Table 4). The hybrid DFT functionals (i.e., B3LYP and LCBLYP) provided closer results. The results from TDHF, which are listed in Table 5, overestimated the excitation energy as compared to the SACCI results. For all the functionals, the DC-TDDFT calculations reproduced the

a

The times required for the preceding HF or DFT calculation are not included. An Intel Xeon X5470 (3.33 GHz) processor was used on a single core.

Table 4. Basis-Set Dependence of the n−π* Excitation Energy (in eV) of C16H17CHO Obtained Using DC- and Conventional CIS, TDDFT, and SACCI Calculationsa

a

basis set

CIS

DC−CIS

(diff)

TDDFT

DC-TDDFT

(diff.)

SACCI

DC-SACCI

(diff.)

6-31G 6-31G** 6-311G 6-311G** cc-pVDZ cc-pVTZ

4.41 4.64 4.41 4.61 4.59 4.63

4.34 4.58 4.31 4.49 4.51 4.58

(−0.07) (−0.06) (−0.10) (−0.12) (−0.08) (−0.06)

3.48 3.53 3.48 3.51 3.49 3.50

3.39 3.50 3.63 3.49 3.46 3.46

(−0.09) (−0.03) (+0.15) (−0.02) (−0.02) (−0.04)

3.85 3.87 3.84 3.79 3.84 −

3.82 3.76 3.80 3.69 3.76 3.78

(−0.03) (−0.11) (−0.04) (−0.10) (−0.08) −

The energy deviations from the conventional results are shown in parentheses. The ES size of nex = 3 was adopted for the DC calculations. 5569



Article

excitation energies obtained using the conventional method within 0.15 eV. 3.3. Photoactive Yellow Protein. The present method was applied to PYP; the crystal structure, which is depicted in Figure 4a, is registered under the accession code 2PHY in the

Table 6. Excitation Energies (in eV) of the HC4 (4hydroxycinnamic acid; pigment), AC (active center; HC4, Tyr42, Glu46, Thr50, Arg52, and Cys69), and entire PYP modelsa CIS TDDFT SACCI Exptl.

HC4

(diff.)

AC

(diff.)

PYP

(diff.)

4.64 4.60 4.52 4.37b

(+0.27) (+0.23) (+0.15)

4.40 3.83 3.59 2.78b

(+1.61) (+1.06) (+0.81)

4.33 3.38 − 2.78b

(+1.55) (+0.60) −

a

The HC4 and AC models were calculated using the conventional method and the PYP model was calculated using the ONIOM scheme. The deviations of the excitation energies from the experimental values are shown in parentheses. bRef 59.

experimental shift; this is due to the lack of interactions with the outside residues. We attempted to include the effects of the outside residues by using the ONIOM method 61 with the Gaussian09 program;45 this method has been widely used for calculating biomolecules in quantum chemistry. The AC model was adopted as the high-level layer where the CIS and TDDFT methods were applied (note that ONIOM-SACCI is not implemented in Gaussian). The other atoms in the PYP system were treated using the universal force field.62 The calculated and experimental excitation energies are shown in Table 6. Although the excitation energies obtained using the ONIOM method are red-shifted from the results of the AC model, they still did not reach the experimental value. The ONIOM method with a classical electrostatic field for the low-layer region also could not reproduce the experimental value. Next, we performed DC-CIS, TDDFT, and SACCI calculations on the C1 and PYP model systems; the results are given in Table 7. The AC region is adopted as the ES.

Figure 4. (a) Structure of PYP. (b) Schematic of the three-model systems: (i) Active center (AC) model [HC4, Tyr42, Glu46, Thr50, Arg52, and Cys69] (blue), (ii) first coordination sphere (C1) model [AC, Ala67, Pro68, Phe96, and Thr98] (green), and (iii) the entire PYP protein (orange).

Protein Data Bank.58 PYP is a photoreceptor protein that comprises 125 amino acid residues and a pigment, which is deprotonated 4-hydroxycinnamic acid (HC4) linked with Cys69 by a thioester bond. PYP shows an absorption maximum at 446 nm (2.78 eV); corresponds to the lowest π−π* excitation energy from HOMO to LUMO, while the absorption maximum of HC4 is at 284 nm (4.37 eV) in water.59 Thus, a large red-shift of 1.59 eV is induced via binding of the pigment to the apoprotein. It is possible to separately quantify the effects of the pigment on the absorption maximum using quantum chemical calculations. We theoretically investigated the red-shift in PYP using DC-CIS, TDDFT, and SACCI calculations. Accordingly, we constructed three models of this system, of which the schematic is depicted in Figure 4b: (i) the active center (AC) model, which comprises HC4, Tyr42, Glu46, Thr50, Arg52, and Cys69, and is depicted in blue; (ii) the first coordination sphere (C1) model, which comprises AC, Ala67, Pro68, Phe96, and Thr98, and is depicted in green; and (iii) the entire PYP protein, which is depicted in orange. The 6-31G basis set was adopted for the following calculations. First, we evaluated the π−π* excitation energy of a HC4 molecule in water solvent. The solvent effect was accounted for using the integral equation formulation of the polarizable continuum model.60 The excitation energies calculated using the conventional CIS, TDDFT, and SACCI methods are given in Table 6 together with the experimental results.59 The differences between the calculated and experimental excitation energies are presented in parentheses. Although the SACCI method gives the best result, with 0.15 eV deviation from the experimental value, the deviations of the CIS and TDDFT results are also small, i.e.,

Novel Approach to Excited-State Calculations of Large Molecules

Recommend Documents