PCM Applied to Ground and Excited State Potential Energy

Article pubs.acs.org/JPCA

DFTB/PCM Applied to Ground and Excited State Potential Energy Surfaces Yoshio Nishimoto* Fukui Institute for Fundamental Chemistry, Kyoto University, 34-4 Takano Nishihiraki-cho, Sakyo-ku, Kyoto 606-8103, Japan S Supporting Information *

ABSTRACT: Accounting for solvent effects in quantum chemical calculations is vital for the accurate description of potential energy surfaces in solution. In this study, we derive a formulation of the analytical first-order geometrical derivative of ground- and excited-state energies within the time-dependent density-functional tight-binding (TD-DFTB) method with the polarizable continuum model (PCM), TD-DFTB/ PCM. The performance of this is then evaluated for a series of halogen-exchange SN2 reactions. DFTB/PCM reproduces DFT results well for isolated monohalogenated methanes, but its agreement for transition structures significantly depends on the halogen element. The performance of TD-DFTB/PCM is evaluated for the excitedstate intramolecular proton transfer (ESIPT) reaction of 3-hydroxyflavone (3HF) in ethanol. TD-DFTB/PCM reproduces the barrier height of the ESIPT reaction in terms of geometry and energy relatively well, but it fails to reproduce the experimental absorption and fluorescence energies as a consequence of the absence of long-range corrections. Computational timings with and without PCM show that the additional cost of PCM for C500H502 is only 10% greater than the corresponding calculation in vacuum. Furthermore, the potential applications of TD-DFTB/PCM are highlighted by applying it to a double-stranded DNA complexed with dye (PDB ID 108D). We conclude that TD-DFTB/PCM single-point calculations and geometry optimizations for systems consisting of more than 1000 and 500 atoms, respectively, is now manageable and that properties predicted with TD-DFTB must be interpreted with care.

■

formulated by Pierotti,20 and repulsion and dispersion free energies, formulated by Floris et al.21,22 Whereas RISM is computationally demanding, another variant of RISM, called RISM-SCF-SEDD,23 can be used to evaluate solvent structures and has recently been used to successfully reproduce trends in excited state properties,24 which PCM could not do. Because the computational costs of ab initio methods become exponentially greater for large systems, efficient semiempirical quantum-mechanical (QM) methods are attracting more attention than ever. One such method, the densityfunctional tight-binding (DFTB) method,25−27 is expected to give reasonable accuracy at a remarkably reduced computational cost. The first generation of DFTB28 is based on the idea of the extended Hückel method, and the self-consistent-charge (SCC) version,29 known as SCC-DFTB or DFTB2, is now used extensively. DFTB2 exploits SCC cycles, which are practically equivalent to SCF cycles, to obtain better accuracy, and makes it possible to take external (e.g., solvent) effects into consideration. Several models have been combined with DFTB: COSMO,30 the generalized Born model,31 the Poisson− Boltzmann,32 and PCM.33 It should be noted that PCM has been combined with other semiempirical methods34,35 and applied to excited states in ref 34 with ZINDO. Despite the

INTRODUCTION An accurate description of solvent effects is indispensable to the rigorous discussion of chemical reactions because a wide variety of these reactions occur in solution. To include solvent effects, one simple and straightforward approach is to explicitly treat solvent molecules, but the number of atoms to be considered in this way is nontrivial. Even with progress in the QM/MM method,1 the computational cost of such calculations remains high. On the other hand, solvent effects can be efficiently considered in an implicit manner through, for example, the generalized Born model,2,3 the polarizable continuum model (PCM),4 the conductor-like screening model (COSMO),5 and the reference interaction site model (RISM).6,7 These first three approaches are known as the dielectric model, and their computational demand is rather modest compared with that of RISM. Among these methods, the PCM approach is widely used and has been combined with a number of computational methodologies: self-consistent-field (SCF) methods (Hartree− Fock and density functional theory (DFT)),8−11 electron correlation methods,12−15 and time-dependent theories.16−19 In PCM, the solute−solvent interface (cavity surface) is discretized into a set of small elements, called tesserae. Each tessera has an apparent surface charge (ASC) that is induced by the electrostatic potential (ESP) of the solute atoms. Solute− solvent electrostatic interactions are calculated with the density of solute and ASCs, and nonelectrostatic interactions are commonly computed as the sum of the cavitation energy, © XXXX American Chemical Society

Received: November 3, 2015 Revised: January 6, 2016

A

DOI: 10.1021/acs.jpca.5b10732 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

the sum of two-center repulsive pairwise contributions, which are independent of electronic structure. The Greek characters μ and ν (later, κ and λ) represent atomic orbitals (AOs). Although the spin index σ (later, τ) is explicitly denoted, only restricted wave functions with integer occupation are used in this study. The zeroth-order Hamiltonian H0μν can be computed with tabulated parameters, whereas the MO coefficients Cμiσ are determined as a solution of the secular equation HC = εSC, where

importance of solvent effects in discussing chemical reactions, it cannot be said that DFTB has, thus far, fully utilized solvent effect tools. In particular, given that only one study of DFTB/PCM33 has been published recently, it is still important to discuss the applicability and performance of DFTB/PCM. In this paper, we first introduce a slightly modified implementation of DFTB/ PCM into the GAMESS-US36 software package. We then, by comparing with DFT calculations, investigate the performance of DFTB/PCM in the calculation of the ground state for a simple SN2 reaction in which the stabilization of halogen ions is considered to be essential to describe potential energy surfaces (PESs) accurately. With regard to DFTB calculations with halogen atoms, DFTB337−39 has recently become more prominent with the publication of a number of reports demonstrating improved parameters.39−43 For this reason, both conventional SCC-DFTB (DFTB2) and DFTB3 will be used and compared in the following sections. As mentioned, to date, solvent effects have been included in implicit ways for DFTB calculations, but there have been very few instances33 of such applications to the time-dependent version of DFTB (TD-DFTB).44 The main aim of this study is to develop the analytical first-order geometrical derivative within TD-DFTB combined with PCM (TD-DFTB/PCM) using the recent implementation of TD-DFTB45 in GAMESSUS. While the previous development of TD-DFTB/PCM by Barone et al.33 omitted the TD-DFTB/PCM gradient, our development includes it and makes it possible to perform efficiently geometry optimizations, including transition structure (TS) searches, and intrinsic reaction coordinate (IRC)46 and vibrational frequency calculations (seminumerical secondorder geometrical derivative) at excited states. TD-DFTB/ PCM is applied to the excited-state intramolecular proton transfer (ESIPT) reaction of 3-hydroxyflavone (3HF) in ethanol. Noting that the advantage of DFTB is its computational efficiency, the computational timings with and without PCM and its application to DNA will also be discussed at the end of this paper.

0 Hμν = Hμν +

qA = =

Eint =

σ

i

μν

ν

μ∈A

(4)

ν

∑ Vkqk̅

(5)

k

where Vk and qk are the ESP exerted by the solute and the induced surface charges, respectively, on tessera k. The latter will always be denoted with a bar on the q, so that we can distinguish the Mulliken charges on an atom and the induced charges on tesserae. The ESP on tessera k, Vk, in this work is defined as

Vk = −∑ A

ΔqA rAk

(6)

where rAk is the distance between atom A and tessera k and Vk is updated in each SCC cycle. The interaction energy can thus be written explicitly within DFTB/PCM as q̅ Eint = −∑ ΔqA ∑ k rAk (7) A k The induced charges qk are obtained by solving the matrix equation

q̅ = −C−1V (8) where C is a matrix that depends on the position of tesserae and the choice of PCM model.4 In this study, C-PCM (conductor PCM)11 is used, and the prefactor of ε/(ε − 1), where ε is the dielectric constant of the solvent, is absorbed into C for simplicity. The ESP on solute atoms, or the contribution to the Hamiltonian, is

where E is the internal energy of the solute, and E is the electrostatic interaction between solute and solvent. The former is defined as 1 2

μ∈A

with the usual density matrix Dμνσ. In eq 1, the electrostatic interaction energy between solute and solvent, Eint, is defined as

(1)

∑ ∑ niσ ∑ CμiσCνiσHμν0 +

i

∑ ∑ ∑ Dμνσ Sμν σ

int

E=

∑ ∑ niσ ∑ ∑ CμiσCνiσSμν σ

METHODOLOGY Overview and Implementation of DFTB/PCM and TDDFTB/PCM. DFTB/PCM has been formulated in ref 33 by Barone et al. in detail, so only a brief explanation is given here. The formulation given here is almost equivalent to the one in ref 33, but there are some differences, such as the way the electron distribution on tesserae is treated, which will be mentioned later. In PCM, the free energy of the solute−solvent system is written as 1 int E 2

(3)

in which Sμν is the overlap matrix element and the last term represents the electrostatic interaction between solute and solvent. The Mulliken charge ΔqA on atom A is the difference between the Mulliken population qA and the number of electrons on neutral atom A. The former is defined as

■

.=E+

1 PCM Sμν ∑ (γAC + γBC)ΔqC + V μν 2 C

∑ ΔqAΔqBγAB + Erep

⎛ 1 1 1 ⎞ PCM = − Sμν ∑ ⎜ + V μν ⎟q 2 r rBk ⎠ k̅ k ⎝ Ak

AB

(2)

where niσ is the electron occupation of the ith molecular orbital (MO). The second term of eq 2 is a Coulomb interaction, and γAB depends on the distance between two atoms, A and B, and their Hubbard values.29 The last term in eq 2, Erep, represents

(9)

which is used for the Hamiltonian (Fock) matrix in eq 3. Approximations introduced in the derivation were discussed in ref 33. Usually, solute charges are decomposed into electrons B


Article


⎛ 1 1 1 ⎞ −1 ⎛ 1 1⎞ + Gμν , κλ = − SμνSκλ ∑ ⎜ ⎟(C )kl ⎜ + ⎟ 4 r rBk ⎠ rDl ⎠ ⎝ rCl kl ⎝ Ak

and nuclei, as are induced charges, but for simplicity, the formulation presented here does not write these contributions separately. The first-order geometrical derivative (gradient) of the energy can be obtained by differentiating the total free energy (eq 1) with respect to the geometrical parameter a, ∂. = ∂a

⎡ ∂H 0

∑ ∑ niσ ∑ CμiσCνiσ ⎢ σ

i

μν

⎢⎣ ∂a

μν

− εiσ

(16)

and that for an MO one is Giaσ , jbτ = −∑

∂Sμν

A

∂a

⎪

⎪

Vkjbτ = −∑

A

k

qAiaσ =

(11)

∑ (A + B)iaσ , jbτ Zjbτ = −R iaσ (20)

jbτ

(12)

with

and

R iaσ =

∑ {(X + Y )ibσ Hab+ σ[(X + Y)]} + Hia+σ[T] b

(13)

−

where elements of the coupling matrix Kiaσ,jbτ are obtained by transforming an AO representation into an MO one. The AO representation based on DFTB2 is 1 SμνSκλ(γAC + γBC + γAD + γBD) 4

∑ {(X + Y )jaσ H ji+σ[(X + Y)]} (21)

j

The vector T (the unrelaxed difference density matrix) is given elsewhere,50 and for an arbitrary vector V′, + H pq σ [V′] =

(14)

∑ (K pqσ ,rsτ + Gpqσ ,rsτ )Vrs′ τ rsτ

for singlet−singlet excitation, and the one for singlet−triplet excitation is Kμν , κλ =

(19)

ν

and are considered to be the surface charges induced by the transition between MOs j and b of τ spin on tessera k. Analytical Gradient for TD-DFTB/PCM. A detailed derivation for TD-DFT/PCM can be found in previous studies by Scalmani et al.18 and by Wang et al.48 Because the derivation of the analytical gradient for TD-DFTB45,49 is somewhat similar to that of TD-DFT/PCM, we avoid repeating it here and simply describe its differences from the analytical gradient without PCM. The analytical gradient within TD-DFTB/PCM requires the solving of Z-vector equations.50,51 To evaluate vector Z, one has to solve the following equation

where matrix elements of A and B are defined as

Kμν , κλ =

iaσ T −1 jbτ

qjbτ k

where μ ∈ A and ν ∈ B, and εiσ is the eigenvalue of the ith MO. ∂r−1 Ak /∂a is the derivative of the inverse distance 1/rAk, and the definition of the term can be found elsewhere.47,48 Note that an actual calculation has to include the displacement of tesserae due to the perturbation along a, and this effect is implicitly included in the derivative term. To obtain excitation energies ω and vectors X and Y, one has to solve the following non-Hermitian equation

Biaσ , jbτ = K iaσ , bjτ + Giaσ , bjτ

(18) jbτ

∑ ∑ CμiσCνaσSμν μ∈A

(10)

Aiaσ , jbτ = δστδabδij(εaσ − εiσ ) + K iaσ , jbτ + Giaσ , jbτ

(17)

so that Giaσ, jbτ = (V ) q = −(V ) C V . Here, it is convenient to define Mulliken atomic transition charges as

−1 ∂rAk 1 ∂C q + q̅ T 2 ∂a ∂a ̅

⎛ 1 0 ⎞⎛ X ⎞ ⎛ A B ⎞⎛ X ⎞ ⎜ ⎟⎜ ⎟ = ω⎜ ⎟⎜ ⎟ ⎝ B A ⎠⎝ Y ⎠ ⎝ 0 −1⎠⎝ Y ⎠

rAk

rAk iaσ T

⎪

∑ ΔqA ∑ qk̅

k

qAjbτ

A

⎪

∂Erep − ∂a

∑

qk̅ jbτ

where qjbτ k are the induced charges obtained by solving eq 8 with the potential on tesserae modified as

1 ∂Sμν ⎧ ⎨∑ (γ + γBC)ΔqC + 2 ∂a ⎩ C AC ⎤ ⎛ 1 ∂γ 1 ⎞ ⎫⎥ 1 − ∑⎜ + ⎟qk̅ ⎬ + ∑ ΔqAΔqB AB rBk ⎠ ⎭⎥⎦ 2 AB ∂a k ⎝ rAk +

qAiaσ

(22)

General MOs are represented with p, q, r, and s indices. Compared to H+pqσ[V′] in a vacuum, the additional term of Gpqσ,rsτ defined in eq 17 has to be evaluated. Note that H−pqσ[V′] is again zero in this study, because the (A−B) matrix is assumed to be diagonal. This is a reasonable assumption given the nature of the PBE functional, which does not contain Hartree−Fock exchange. Very recent progress in DFTB can include long-range corrections52,53 with which H−pqσ[V′] must be evaluated; however, the DFTB model used in the present study does not include them. After solving the Z-vector equation, we construct the relaxed one-particle difference density matrix (P = T + Z) and the Lagrange multiplier (W; the energy-weighted difference density matrix), which are defined in ref 50. With these quantities, the excitation energy gradient is

1 SμνSκλ(MAδAC + MBδBC + MAδAD + MBδBD) 4 (15)

where μ ∈ A, ν ∈ B, κ ∈ C, and λ ∈ D, and MA represents the spin constant of atom A, which is also obtained from reference DFT calculations. Note that i and j indicate occupied orbitals, while a and b indicate virtual orbitals. Although the formulation given here is for DFTB2, the extension to DFTB337−39 requires no additional changes in the connection with PCM. See ref 45 for the MO representation and the details of TD-DFTB3. The PCM contributions are expressed with the Giaσ,jbτ term in eqs 12 and 13. The expression for an AO representation is C


Article

The Journal of Physical Chemistry A ∂ω = ∂a

0 ∂Hμν

∑

∂a

μνσ

∂a

μνσκλτ vac

=

∂ω ∂a

+

∑

∂Sμν

μνσ

∂Kμν , κλ

∑

+

Pμνσ −

∂a

∂Gμν , κλ

∑

̂ , κλτ + Γμνσ

∂a

μνσκλτ

∂Gμν , κλ

∑

∂a

μνσκλτ

where Rk and RA are coordinates of tessera k and atom A, in ref 33. The different definition comes from the fact that the previous DFTB/PCM implementation employed a Gaussian charge distribution for tesserae to achieve a smooth solvation potential.58 The proposed scheme expands surface charges on tesserae with a Gaussian function

Wμνσ ̂ , κλτ Γμνσ

̂ , κλτ Γμνσ

⎛ ζ 2 ⎞3/2 ϕk (r;R k) = ⎜ k ⎟ exp( −ζk 2|r − R k|2 ) ⎝ π ⎠

(23)

where ∂ωvac/∂a is the gradient without explicit PCM contributions, and its lengthy expression within TD-DFTB is given in ref 45. The effective two-particle difference density matrix Γ̂μνσ,κλτ is ̂ , κλτ = (X + Y )μνσ (X + Y )κλτ + PμνσDκλτ Γμνσ

where the exponent ζk is determined from the angular quadrature weights and a constant.58 After mathematical manipulation, it can be shown that

(24)

⟨A|k⟩ =

50

as usual without the Hartree−Fock exchange. The gradient in eq 23 can be simplified by following the procedure in refs 18 and 48 using the notation for qPA and qXY A in ref 45. The working expression is ∂ω ∂ωvac 1 = − ∂a ∂a 2 −

−

1 2

μν

∂Sμν

A

∂a

∂a

− 2 ∑ qAXY

∑

A

k

⎛ 1 1 ⎞ P + ⎟q rBk ⎠ k̅ ⎝ rAk

k

∂Sμν ∂a

∂r −1 qk̅ P Ak

k

⎛ 1 1 ⎞ + ⎟q rBk ⎠ k̅ ⎝ rAk

∑⎜

μν

∑ ΔqA ∑

k

−

⎛ 1 1 ⎞ XY + ⎟q rBk ⎠ k̅ ⎝ rAk

∑⎜ k

∑ qAP ∑ qk̅ A

k

−1 ∂rAk ∂a

∂r −1 qk̅ XY Ak ∂a

∂C ∂C XY + (q̅ P)T q̅ + (q̅ XY )T q ∂a ∂a ̅

(25)

In eq 25, qPk and qXY k are the ASCs induced by P and the excitation vector (X + Y), respectively. Implementation of (TD-)DFTB/PCM. TD-DFTB/PCM was implemented in GAMESS-US,36 which contains already developed DFTB45,54,55 and PCM48,56,57 modules. All equations presented in previous subsections are based on DFTB2,29 but the extension to DFTB3 requires only the modification of the internal energy, with the solute−solvent interactions remaining the same. This is because the ESP on the tesserae are computed through Mulliken charges on solute atoms with eq 6, which are indirectly affected by the choice of the DFTB model through the change of electronic structure. PCM contributions to the Hamiltonian (eq 9) are decided by geometrical parameters and induced charges. Some quantities defined in this work are different from those in the work by Barone et al.33 For instance, eq 7 in the present work should be read as (see eqs 8, 25, and 26 in ref 33) Eint = −∑ ΔqA ∑ ⟨A|k⟩qk̅ A

(26)

k

with ⟨A|k⟩ =

∬

ϕk (r;R k) δ(r′ − RA) |r − r′|

dr dr′

erf(ζk|RA − R k|) |RA − R k|

(29)

Our definition, ⟨A|k⟩ = 1/|RA− Rk| = 1/rAk, is thus missing the error function of the numerator. Nevertheless, we use the simpler definition (⟨A|k⟩ = 1/rAk) in our implementation. In GAMESS-US, the electrostatic interactions between tesserae are described with a simple-point charge model. Considering that the integral in eq 29 is the consequence of Gaussian functions on each tessera, eq 29 should not be used in GAMESS-US unless the PCM part is modified extensively. To be consistent with the PCM implementation in GAMESS-US, ⟨A|k⟩ = 1/rAk should be the definition used. Note that there are some other differences in the implementation of PCM in GAMESS-US. For instance, by default, GAMESS-US employs the fixed points with variables areas (FIXPVA)57 tessellation scheme in defining the cavity surface. One can choose other schemes to generate tesserae, but FIXPVA is considered to be the best among those implemented in GAMESS-US in terms of computational cost and the accuracy of the gradient. As a comparison with the previous implementation in ref 33, the solvation energies for 18 neutral small molecules (see Table 2 in ref 33; α = 1.1) are compared. The root-mean-square (RMS) and maximum deviations in solvation energies in ref 33 and our implementation are 0.25 and 0.68 kcal/mol with DFTB2. A similar comparison with PBE/STO-3G gives deviations of 0.33 and 0.55 kcal/mol. Thus, these deviations may be largely attributed to differences in the implementation of PCM. In addition, excitation energies of Nile Red in a vacuum, heptane, benzene, and acetonitrile show that values calculated with GAMESS-US are higher than those in ref 33 by 0.04 eV for all entries. As the excitation energy in a vacuum is slightly overestimated with our code, the systematic overestimation should again be attributed to differences in the implementation. Note that the DFTB used in ref 33 uses analytical expressions, called “DFTBA”,59 while DFTB in GAMESS-US uses tabulated parameter files. To verify the accuracy and the correctness of the implemented analytical gradient, we compared the analytical and numerical gradients at the ground state and the three lowest singlet−singlet and singlet−triplet excited states of C2H5NO2 (glycine) produced by ChemCraft.60 As shown in Table S1 (see the Supporting Information), the error of gradient is less than 1.0 × 10−8 with both DFTB2/PCM and DFTB3/PCM at the ground state. The largest deviation, of 9.7 × 10−8, is found for the first singlet−singlet excitation with TDDFTB3/PCM, which is still fairly small.

∑⎜

∂a

μν

∑ ∑ (X + Y )μνσ σ

−

σ

∑ ∑ Dμνσ σ

∂Sμν

∑ ∑ Pμνσ

(28)

(27) D



■

Article

■

COMPUTATIONAL DETAILS We used the MIO29,39,61(+halorg62) and 3OB40,41,43 sets of DFTB parameters63 in DFTB2 and DFTB3 calculations, respectively. Hubbard derivatives and the ζ parameter (4.0) of the damping function for hydrogen atoms in DFTB3 calculations were also taken from the references above. A tighter convergence criterion for SCF (CONV = 1.0D-08 in $SCF) was used in all computations. Geometry optimizations were performed until the RMS and maximum gradient elements became smaller than 1/3 × 10−6 and 10−6 Hartree/ Bohr, respectively (OPTTOL = 1.0D-06 in $STATPT). The FIXPVA57 tessellation scheme was used, and ASCs, solutions of eq 8, were obtained with an iterative solver64 instead of the direct inversion of C. To minimize the effect of rotational variance, the number of initial tesserae per atom was increased from 60 to 240. Although GAMESS-US can calculate nonelectrostatic contributions (cavitation, repulsion, and dispersion free energies) for DFTB, they were not included in any of the calculations presented hereafter. Second-order geometrical derivatives were numerically computed with analytical first-order geometrical derivatives. In DFT calculations, a finer Lebedev grid was used (NRAD = 99 and NLEB = 590 in $DFT). DFTB/PCM was applied to SN2 reactions in water. DFTB3 calculations included the D3(BJ) correction65,66 with suitable parameters43 for systems that contain halogens. No dispersion corrections were added to DFTB2 calculations, because they are not considered in the parametrization.62 In fact, the effect of the dispersion corrections using parameters optimized for DFTB3 is not significant, as shown in Table S2 (see the Supporting Information), and the difference in relative energy is less than 1 kcal/mol. As a reference, B3LYP67,68-D3(BJ) and PBE69,70-D3(BJ) calculations using the aug-cc-pVTZ(-PP)71−75 and SBKJC(d)76,77 basis sets were performed. The aug-ccpVTZ-PP functions for Br and I were taken from EMSL Basis Set Exchange.78,79 A pseudopotential (denoted by a suffix of “-PP”)74,75 was utilized for Br and I, and 10 and 28 electrons were assigned as core electrons, respectively. SBKJC(d) indicates that the SBKJC basis set for halogen atoms was polarized with a d-polarization function (exponent = 0.8, coefficient = 1.0), while the original SBKJC was used for C and H. The SBKJC(d) basis set for halogen was a [2s,2p,1d] contracted basis with a pseudopotential, which is close to the quality of the “basis set” used in DFTB calculations. The accuracy of DFTB at excited states is discussed with 3hydroxyflavone (3HF). The PES of the ground and the first excited state and the absorption and fluorescence energies in ethanol are discussed with (TD-)DFT and (TD-)DFTB results. Geometries were optimized with (TD-)DFTB2, (TD-)DFTB3, and (TD-)B3LYP/6-31G(d,p).80,81 Energy refinement with B3LYP, PBE, BLYP,82,83 and LC-BLYP (μ = 0.33)84 functionals using the aug-cc-pVDZ basis set71,72 was performed for optimized geometries at (TD-)B3LYP/6-31G(d,p). Dispersion corrections were not included for 3HF. As a demonstrative example, (TD-)DFTB3/PCM was applied to the fluorescence enhancement of a dye molecule upon binding to DNA. The Slater−Kirkwood dispersion correction85 was used, as the parameters for this dispersion correction have been optimized for DNA. The threshold of geometry convergence and the number of initial tesserae were set to the default values (OPTTOL = 1.0D-04 and NTSALL = 60, respectively).

RESULTS AND DISCUSSION Accuracy of DFTB: SN2 Reaction. Let us briefly discuss an elemental reaction in chemistry: the SN2 reaction. This involves an ionic fragment, meaning that solvent effects are important in accurately describing its PESs. In this study, a set of halogen exchange reactions (Figure 1) in water was considered.

Figure 1. A representative SN2 reaction. X1 = X2 = F, Cl, Br, and I in this study.

Although higher level calculations, such as CCSD(T) with PCM, would be desirable as a reference, such calculations are not possible with GAMESS-US. Therefore, a widely used level of theory served as a reference instead: B3LYP-D3(BJ)/aug-ccpVTZ(-PP). Table 1 shows the C−X distance (RC−X; X = F, Cl, Table 1. Calculated C−X (X = F, Cl, Br, and I) Distances (RC−X in Å) for the Transition Structures (TSs) and Isolated CH3X and ΔG (in kcal/mol) in Water with PCM RC−X method DFTB2 DFTB3-D3(BJ) PBE-D3(BJ) B3LYP-D3(BJ) PBE-D3(BJ) B3LYP-D3(BJ) DFTB2 DFTB3-D3(BJ) PBE-D3(BJ) B3LYP-D3(BJ) PBE-D3(BJ) B3LYP-D3(BJ) DFTB2 DFTB3-D3(BJ) PBE-D3(BJ) B3LYP-D3(BJ) PBE-D3(BJ) B3LYP-D3(BJ) DFTB2 DFTB3-D3(BJ) PBE-D3(BJ) B3LYP-D3(BJ) PBE-D3(BJ) B3LYP-D3(BJ)

parameter/basis set X=F MIO 3OB SBKJC(d) SBKJC(d) aug-cc-pVTZ aug-cc-pVTZ X = Cl MIO 3OB SBKJC(d) SBKJC(d) aug-cc-pVTZ aug-cc-pVTZ X = Br MIO 3OB SBKJC(d) SBKJC(d) aug-cc-pVTZ-PP aug-cc-pVTZ-PP X=I MIO 3OB SBKJC(d) SBKJC(d) aug-cc-pVTZ-PP aug-cc-pVTZ-PP

for CH3X

for TS

ΔG at TS

1.407 1.383 1.436 1.437 1.414 1.409

1.930 1.924 1.880 1.882 1.872 1.875

26.22 31.09 17.50 21.37 24.91 28.35

1.820 1.804 1.834 1.841 1.802 1.808

2.275 2.355 2.376 2.398 2.329 2.353

24.18 22.71 17.75 19.34 20.60 22.83

1.979 1.970 1.997 2.003 1.966 1.972

2.420 2.593 2.524 2.544 2.476 2.498

21.04 31.47 15.28 17.02 16.70 19.17

2.183 2.166 2.196 2.200 2.165 2.171

2.625 2.662 2.730 2.750 2.674 2.697

19.02 16.85 13.45 15.32 14.32 17.00

Br, and I) for isolated CH3X and transition structures (TSs) and ΔG at TSs. The reference of ΔG was set to the sum of G for isolated X− and CH3X, because there are no ion−dipole complexes (see, for instance, refs 86, 87, and 88) in solution phase. DFTB3 parameters involving halogen atoms were parametrized43 to reproduce minimum geometries optimized with B3LYP/def2-TZVPP, so we can expect DFTB3-D3(BJ) to E


Article


Accuracy of DFTB: Excited-State Intramolecular Proton Transfer. 3HF is known to exhibit an excited-state intramolecular proton transfer (ESIPT) reaction.89 Because of the simplicity and the availability of experimental90 absorption and emission properties, the ESIPT reaction of 3HF should be suitable for testing (TD-)DFTB/PCM here. In the manner described in Figure 2, the enol form of 3HF is dominant at the

reproduce the RC−X obtained with B3LYP-D3(BJ)/aug-ccpVTZ-(PP) well. This is indeed so, according to the RC−X for isolated CH3X. The deviation of RC−X for X = Cl, Br, and I is less than 0.01 Å, and the largest deviation, for X = F, is only 0.026 Å. These deviations are comparable with PBE-D3(BJ)/ aug-cc-pVTZ(-PP) results, except for X = F, and they are smaller than those obtained with the SBKJC(d) basis set. As for the distance for X = F, DFTB2 somehow reproduces the B3LYP result very well, but DFTB3-D3(BJ) does better than DFTB2 for other halogen elements. Table S3 (see the Supporting Information) compares the bond lengths in vacuum and in water. RC−X in water is longer than that in vacuum by 0.02 Å for X = F and by 0.01 Å for X = Cl, Br, and I, indicating that solvent effects have a small effect on RC−X. Considering that reference molecules in the parametrization are PES minima, the applicability of DFTB with parameters generated in this way for TSs of the SN2 reaction is not clear. DFTB3-D3(BJ) does a relatively good job of reproducing RC−X and ΔG obtained with B3LYP-D3(BJ)/aug-cc-pVTZ(-PP) for X = Cl and I, and it is actually better than the DFT calculations done with the SBKJC(d) basis set. RC−X for X = Cl obtained with DFTB3-D3(BJ) deviates by only 0.002 Å from B3LYPD3(BJ)/aug-cc-pVTZ, and for X = I it is acceptably small, 0.035 Å. DFTB3-D3(BJ) reproduces ΔG obtained with B3LYPD3(BJ)/aug-cc-pVTZ(-PP) within 0.2 kcal/mol for X = Cl and I and, thus, satisfactorily succeeds in reproducing B3LYP results. However, DFTB3-D3(BJ) has a deviation of 0.049 Å for X = F and 0.096 Å for X = Br in RC−X. DFTB2 reproduces RC−Br slightly better, with a deviation of 0.077 Å, but this is still large, and DFTB2 performs worse for other halogen elements. We tried to identify the reason why DFTB3-D(BJ) gives a larger discrepancy for the RC−Br distance than DFTB2, by checking, one by one, the following differences between the two methods: (1) the third-order term of a Taylor expansion of the exchange−correlation energy, (2) the damping function for hydrogen atoms, (3) the parameter set (MIO versus 3OB), and (4) the dispersion correction [D3(BJ)]. Of these, the addition of differences 1, 2, and 4 did not have a significant influence on RC−Br, which indicates that the difference in the parameter set is the dominant factor for the discrepancy of RC−Br. In fact, RC−Br obtained with DFTB3-D3(BJ) using the MIO parameter set is 2.419 Å, which is very close to the DFTB2 result. Although mixing two parameter sets is physically meaningless, if we replace the MIO parameters for C−Br and Br−C with those of 3OB, RC−Br for the TS becomes 2.506 Å. Considering that PCM has not been used in the parametrization, we compared RC−Br obtained in a vacuum; DFTB3-D3(BJ) again overestimated RC−Br by 0.096 Å. It appears that the overestimation of RC−Br for the TS is mostly a result of the 3OB parameters, rather than the DFTB3 model or dispersion effects. Regarding energetics, the ΔG of DFTB for X = F deviates by a few kcal/ mol compared with B3LYP-D3(BJ)/aug-cc-pVTZ-PP, but DFTB3-D3(BJ) gives a significant deviation of more than 10 kcal/mol for X = Br. Even if we use the MIO parameter set, other than for C−Br and Br−C in the DFTB3-D3(BJ) calculation, as mentioned above, ΔG becomes worse: 35.86 kcal/mol. Considering these facts, the MIO set seems to be better than the 3OB set only for X = Br in the SN2 reaction. As a conclusion, the applicability of DFTB3 to the prediction of TSs of the halogen exchange reaction strongly depends on the halogen elements involved. DFTB3-D3(BJ) performs well for Cl and I, but care is needed for F and Br.

Figure 2. Tautomerization of 3HF initiated by light absorption. The wavelengths are taken from experiment.90

ground state (S0). Upon absorption of light, the proton on Oa in the enol form is easily transferred onto Ob through the ESIPT reaction at the first singlet−singlet excited state (S1), and the keto form becomes energetically more favorable than the enol form at S1. 3HF has dual emission character caused by an emission at the Franck−Condon region and the keto form with short and long wavelengths (402 and 533 nm in ethanol90), respectively. The energies of the TS and keto form with TD-DFTB and TD-DFT calculations at S0 and S1, relative to the enol form, are summarized in Table 2. In contrast to the results for the SN2 reaction discussed in the previous subsection, both DFTB2 and DFTB3 reproduce the TS energy obtained with B3LYP/6-31G(d,p) at both S0 and S1 within a few kcal/mol deviation in terms of ΔE. Interestingly, ΔE for the TS and keto form calculated with B3LYP/6Table 2. Calculated Relative Energies (ΔE in kcal/mol) of 3HF at the Ground (S0) and the First Excited (S1) States with Respect to Enol Formsa method DFTB2 DFTB3 B3LYP/6-31G(d,p) B3LYP/aug-cc-pVDZb PBE/aug-cc-pVDZb BLYP/aug-cc-pVDZb LC-BLYP/aug-cc-pVDZb TD-DFTB2 TD-DFTB3 TD-B3LYP/6-31G(d,p) TD-B3LYP/aug-cc-pVDZb TD-PBE/aug-cc-pVDZb TD-BLYP/aug-cc-pVDZb TD-LC-BLYP/aug-cc-pVDZb

TS

keto

S0 14.24 (12.11) 10.71 (8.63) 12.17 (10.44) 12.65 (10.93) 8.16 (6.43) 10.09 (8.36) 14.08 (10.61) S1 2.12 (0.23) 2.08 (0.33) 1.79 (0.07) 3.27 (1.54) 1.21 (−0.52) 2.51 (0.79) 1.47 (−0.25)

12.32 (11.50) 9.14 (9.47) 10.91 (10.91) 10.44 (10.44) 7.27 (7.56) 7.97 (8.26) 12.42 (12.87) −10.71 (−9.24) −7.21 (−7.03) −9.66 (−9.16) −8.24 (−7.75) −6.89 (−6.39) −7.27 (−6.78) −8.81 (−8.31)

a Values in parentheses indicate relative free energies (ΔG in kcal/ mol). bGeometries were optimized at (TD-)B3LYP/6-31G(d,p).

F


Article


Figure 3. Relative energies (in kcal/mol) with respect to the enol form during IRC scan projected onto the Oa−H distance (see Figure 2) in Å at (A) S0 and (B) S1. Black circle, blue tilted square, and red square symbols indicate (TD-)DFTB2, (TD-)DFTB3, and (TD-)B3LYP/6-31G(d,p) results, respectively, and blank symbols correspond to TSs.

31G(d,p) at both S0 and S1 lies between the DFTB2 and DFTB3 energies. Considering that the relative energies obtained with four different functionals (B3LYP, PBE, BLYP, and LC-BLYP) have a spread of more than 2 kcal/mol and those obtained with DFTB3 are within the error bar at both S0 and S1, it is reasonable to conclude that DFTB3 gives a good description of the PES for the ESIPT reaction, and is certainly comparable with DFT. The DFTB2 energetics for S0 is very close to the LC-BLYP/aug-cc-pVDZ//B3LYP/6-31G(d,p) results. Figure 3 illustrates the PESs at S0 and S1 with (TD-)DFTB2, (TD-)DFTB3, and (TD-)B3LYP/6-31G(d,p). The energies were obtained with the intrinsic reaction coordinate (IRC)46,91 calculation, and the geometries along the IRC paths were projected onto the Oa−H bond distance (see Figure 2). At S0, the B3LYP result is between those of DFTB2 and DFTB3, but DFTB3 reproduces the B3LYP result slightly better, as can be seen from the values given in Table 2. The Oa−H distances at stationary points, which are shown in Table S4 (see the Supporting Information), also support this observation. From a comparison of the barrier heights (Table 2) and the Oa−H distances (Table S4) at the TS it can be seen that the S1 PESs around the TS for the three methods are seemingly similar. However, for the IRC calculation with TD-DFTB2, the S1 PES crosses (Figure 4A) the S2 PES at around 1.78 Å. The S1 state is characterized as a transition from the highest occupied MO (HOMO) to the lowest unoccupied MO (LUMO) around the TS, and the S2 is a transition from the HOMO−1 (second highest MO) to the LUMO (see Figure 4B). In contrast, for the IRC toward the keto form, the frontier orbitals (Figure 4C) are similar to those in the TS, but the S2 PES is lower than the S1 PES. The oscillator strength of the HOMO−1 → LUMO excitation is close to zero, so the crossing between S1 and S2 is likely a drawback of (TD)DFTB2, considering that 3HF has dual-emission character. All the TD-DFTB2 results were obtained by choosing the HOMO → LUMO excitation during the IRC scan. The PESs obtained with TD-DFTB3 and TD-B3LYP do not have such a crossing, so (TD-)DFTB3 seems to be a better choice for the present example. Absorption energies can be computed using the linearresponse nonequilibrium approach17 with GAMESS-US. Note that the corrected linear-response92 or state specific model93,94 has not been implemented in GAMESS-US. As these sophisticated approaches are not currently available, fluorescence energies are simply obtained as a difference between free energies (eq 1) at equilibrated S1 and S0, with a static dielectric

Figure 4. (A) S1 and S2 PESs during IRC scan projected onto Oa−H distance with TD-DFTB2. HOMO−1, HOMO, and LUMO (±0.03 (e/a0)1/2 isosurface) of the (B) TS and (C) keto form. Values in parentheses are the eigenvalues of the vector in Hartree.

constant of 24.55. This is a crude approximation, but there is no need to obtain the excited state density once or iteratively. Table 3 summarizes calculated and experimental90 absorption and fluorescence energies. Since 3HF in ethanol has a dualemission character, Table 3 has two types of fluorescence energies that correspond to the emissions at the Franck− Condon region (enol) and the tautomerized (keto) forms. Unfortunately, TD-DFTB/PCM underestimates the experimental absorption and fluorescence energies by more than 0.4 eV. Considering that generalized gradient approximation (GGA) functionals (PBE and BLYP) also predict lower excitation energies than hybrid (B3LYP) or long-range corrected (LC-BLYP) functionals, it is reasonable that the significant deviation of PBE, BLYP, and DFTB is largely attributed to the absence of the (long-range) Hartree−Fock exchange. Interestingly, both DFTB2 and DFTB3 reproduce the absorption energies obtained with PBE and BLYP functionals within a deviation of 0.05 eV, and all these methods G


Article

The Journal of Physical Chemistry A Table 3. Calculated and Experimental90 Absorption and Fluorescence Energies in eVa

additional cost is only 10% of the computational cost of the TD-DFTB in a vacuum. Further investigation proves that the most time-consuming step in the calculations with PCM is the iterative computation of ASCs,64 specifically the potential between tesserae (named ASCPOT in GAMESS-US). For example, the sum of the timings for ASCPOT during GS for C500H502 is 25.7 s and corresponds to 91% of the additional timing due to PCM. Computations of ESP between solute and solvent (eqs 6 and 9) are relatively fast. Nevertheless, the iterative solver64 for large systems is considered to be much more efficient than the direct inversion, since the latter requires directly obtaining the C−1 matrix with dimensions of NTS × NTS, where NTS is the number of tesserae, and NTS is 14 070 for C500H502. As polyacetylene is nonpolar, iterative ASC equations (eq 8) for polyacetylene are expected to quickly converge. We performed a single-point calculation for the hen egg-white lysozyme (PDB ID 1IO5)95 and compared the performances. 1IO5 contains a zwitterionic terminal and 18 positively and 9 negatively charged residues, giving a total charge of +9. Although the 1IO5 protein is rather large (1961 atoms with 5014 basis functions), NTS is 11 806, which is comparable with that of C400H402 (NTS = 11 270). The SCF convergence criteria were satisfied after 22 (1IO5) and 17 (C400H402) SCF iterations, whereas 483 and 320 iterations were required to obtain ASCs in total by solving eq 8. If the latter numbers are normalized with the numbers of SCF iterations, the normalized ASC iteration for 1IO5 (22.0) is bigger than that for C400H402 (18.8). The convergence of the ASC equations for polar systems is thus indeed more demanding than that of nonpolar systems. Application of TD-DFTB/PCM. The efficiency of the TDDFTB/PCM method should allow for larger systems to be tackled. To demonstrate its use, we applied the method to a complex of double-stranded DNA with a bis-intercalating dye, namely, d(CGCTAGCG)2 with TOTO dye (PDB ID 108D).96 Experiment96−99 has shown that the intensity of fluorescence of TOTO is greatly enhanced (>3000-fold) upon binding (intercalating) to d(CGCTAGCG)2. Figure 5A shows the optimized ground-state structure of the complex with DFTB3/ PCM. The TOTO, shown with spheres, is bis-intercalated in the minor groove of the d(CGCTAGCG)2, shown with lines. 108D consists of 619 atoms (504 and 115 atoms for DNA and TOTO, respectively), so direct computation with conventional TD-DFT or other wave-function methods is extremely challenging. In addition, without solvent effects, we could not achieve SCF convergence with DFTB. The total charge of TOTO is +4 (Figure 5B), so the 108D complex has an overall charge of −10.

fluorescence method DFTB2 DFTB3 B3LYP/6-31G(d,p) B3LYP/aug-cc-pVDZb PBE/aug-cc-pVDZb BLYP/aug-cc-pVDZb LC-BLYP/aug-cc-pVDZb experiment a b

absorption 3.02 3.04 3.52 3.44 3.07 3.06 3.86 3.59

(410) (408) (352) (360) (404) (405) (321) (345)

enol 2.46 2.52 2.97 2.89 2.62 2.61 3.19 3.08

keto

(505) (493) (418) (430) (473) (474) (389) (402)

1.46 1.90 2.17 2.20 2.04 2.03 2.34 2.33

(850) (654) (571) (564) (609) (611) (530) (533)

Values in parentheses are the corresponding wavelengths (in nm). Geometries were optimized at B3LYP/6-31G(d,p).

underestimate the experimental absorption energy by more than 0.5 eV. The prediction of the absorption energy with DFTB is comparable to that of the GGA functionals. The same is true for the fluorescence energies: DFTB3 underestimates the experimental values by more than 0.4 eV, and the GGA functionals do so by 0.3 eV, too. In spite of an overestimation of the absorption energy with a long-range corrected functional, LC-BLYP, it does give a good prediction of the fluorescence energies. Therefore, long-range corrections clearly improve the prediction of the fluorescence energies. This fact strongly suggests that long-range corrections for DFTB52,53 are necessary to obtain good accuracy with TDDFTB. Computational Timings. Since PCM requires additional computation and how expensive it is was not discussed in ref 33, the comparison of computational timings with and without PCM is worth mentioning. As a test case, we chose a series of trans-polyacetylene (C100nH100n+2; n = 1−5), which contained 500n + 2 basis functions. The five lowest excitation energies were calculated, and the first excitation vectors were used in the subsequent gradient calculations. The initial number of tesserae per atom was decreased to a default value (60). Table 4 shows the CPU timings of TD-DFTB3 calculations with and without PCM (water). Total CPU timings are divided into three major steps: ground state energy (GS), excitation energy (EE), and gradient calculation (GRAD). It can be seen that the inclusion of PCM contributions is a heavy burden for smaller systems. The ratios of the PCM to the vacuum timings are 2.42, 1.80, 1.44, 1.31, and 1.10 for n = 1, 2, 3, 4, and 5, respectively. This means that the additional cost due to PCM is 1.42 times greater than that of the TD-DFTB calculation in a vacuum for the smallest system consisting of 202 atoms. In contrast, for the biggest system, consisting of 1002 atoms, the

Table 4. CPU Timings (in s) To Obtain the Ground State Energy (GS), Five Lowest Excitation Energies (EE), and the Gradient of Energy (GRAD) at the First Excited State for C100nH100n+2 (n = 1−5) with TD-DFTB3/PCM Using the 3ob-3-1 Parameter Seta step

n=1

n=2

n=3

n=4

n=5

GS EE GRAD miscellaneous total

3.8 (2.5) 12.8 (5.0) 2.8 (0.5) 0.2 (0.1) 19.6 (8.1)

20.7 (15.4) 83.2 (44.0) 10.3 (4.2) 0.3 (0.1) 114.4 (63.7)

76.2 (64.7) 255.7 (169.9) 27.2 (14.4) 0.7 (0.4) 359.8 (249.4)

206.6 (187.4) 656.2 (476.3) 63.9 (42.5) 1.5 (1.1) 928.2 (707.3)

444.3 (416.2) 1645.3 (1490.1) 139.2 (103.7) 4.3 (3.0) 2233.1 (2013.0)

a

Calculations were performed with one CPU core of a Xeon E5-1650 v3 (3.50 GHz), and the values in parentheses are CPU timings without PCM (in a vacuum). H


Article


Figure 5. (A) Complex of d(CGCTACG)2 (shown by the lines) with TOTO (spheres) and (B) chemical formula of TOTO.

We first discuss excitations of the free (unbound) TOTO. According to TD-DFTB3/PCM calculations, the absorption should occur for the transition to S3 (HOMO → LUMO) and S4 (HOMO−1 → LUMO+1) states. In contrast, S1 (HOMO → LUMO+1) and S2 (HOMO−1 → LUMO) are dark states with zero oscillator strength. These orbitals are presented in Figure S1 (see the Supporting Information), and the benzothiazole and quinoline groups mainly absorb light for the S3 and S4 excitations. Our TD-DFTB3/PCM calculation predicts that the absorption maxima for TOTO lie at 477.2 and 471.3 nm, which are close to the experimental maximum, 481 nm,97 for free TOTO. We used U-shaped TOTO, but the same calculation with straight TOTO gives almost the same result, with only a trivial difference. We then performed geometry optimization on the S3 and S4 PESs. In both cases, we observed rotation around the double bond between quinoline and the carbon between benzothiazole and quinoline groups. As a result, the relative orientation of the benzothiazole group became perpendicular to the quinoline, and the excitation character turns into a charge-transfer-type excitation, the oscillator strength of which is only 1.8 × 10−5. Because TOTO dye contains 105 atoms with two sets of benzothiazole and quinoline groups, further analysis is somewhat complicated. Instead, we analyzed the character of a TO molecule (Figure 6A) in detail. Note that as indicated in Figure 6A, there are resonance structures; the benzothiazole and quinoline groups are actually conjugated. The molecule has one set of benzothiazole and quinoline groups, so it should be a good model of TOTO. According to the observation above with TD-DFTB3/PCM, the S1 PES is closely related to the twist around the double bond (ψ in Figure 6A, left). We first discuss calculations with the CAM-B3LYP100 functional. Geometry optimization with (TD-)CAM-B3LYP/6-31G(d,p) found one minimum on the S0 PES with ψ = −16.2° and two minima on the S1 PES with ψ = −40.0° and −90.6°. With CAM-B3LYP/aug-cc-pVDZ single-point calculations for the optimized geometry, the absorption maximum is at 425.4 nm and the emission maxima are at 520.9 (ψ = −40.0°) and 678.6 (−90.6°) nm. The oscillator strength for the absorption is 1.014, while those for emission are 0.897 (ψ = −40.0°) and 4.40 × 10−3 (ψ = −90.6°), respectively. Considering that weak emission of free TO is observed experimentally, the lifetime of the minimum at ψ = −40.0° is short, and the population of the structure with ψ = −90.6° should be dominant. Although the absorption maximum for TO is somehow underestimated,

Figure 6. (A) Structure of TO and the definition of the dihedral angle ψ. Its frontier orbitals [±0.04 (e/a0)1/2 isosurface] for (B) S0 and (C) S1 minima with DFTB3/PCM. For parts B and C, values in parentheses represent the eigenvalues of vectors in Hartree.

compared to the experimental value (501 nm97), the experimental emission maximum (640 nm97) is better reproduced with ψ = −90.6°. It seems that this reduction of emission is attributed to the rotation around ψ on the S1 PES. The relative energies of the two minima on the S1 PES with CAM-B3LYP/aug-cc-pVDZ are 59.34 (ψ = −40.0°) and 59.29 (−90.6°) kcal/mol with respect to the S0 minimum. The Franck−Condon structure (ψ = −16.2°) is less stable than these minima with a relative energy of 67.21 kcal/mol at S1. There must be (at least) one TS between the two minima, and the estimated barrier height on the S1 PES through a relaxed scan around ψ is 3−4 kcal/mol with ψ ≈ − 50.0°, so this small barrier can be easily overcome. The twisted structure with ψ = −90.6° is thus accessible on the S1 PES, and the comparison of the emission maximum for experiment and calculation implies that the population of the twisted structure is expected to be dominant. Note that the lambda diagnostic value for the twisted structure is 0.211; therefore, the excitation is characterized as highly charge transfer dependent. (TD-)DFTB3/PCM identified one minimum on both S0 and S1 with ψ = −3.0° and −89.8°. The DFTB calculation thus could not identify a half-rotated minimum on S1. The relative energies of the Franck−Condon structure and the S1 minimum with respect to the S0 minimum are 60.58 and 37.95 kcal/mol, so the twisted structure is fairly favorable and spontaneously found on the S1 PES. TD-DFTB3/PCM predicts that the absorption and emission maxima lie at 472.0 and 1419.4 nm, respectively, whereas the experimental values are 501 and 640 nm.97 The experimental absorption maximum is reasonably reproduced with TD-DFTB3/PCM, but the emission is greatly overestimated, or the emission energy is underestimated by 1.06 eV. The discrepancy of the emission maximum with DFTB is attributed to the absence of long-range corrections. As the CAM-B3LYP calculation indicates, the emission for the twisted structure involves charge transfer, and it is noteworthy that a recent study53 showed that long-range corrections can reduce error for charge transfer excitations by more than 2 eV. The oscillator strengths for the absorption and emission are 0.886 and 6.4 × 10−5, so DFTB predicts that the rotation around ψ gives red-shifted emission with a smaller oscillator strength, in I


Article

The Journal of Physical Chemistry A qualitative agreement with the observation from the CAMB3LYP calculations. The frontier orbitals obtained with DFTB3/PCM (Figure 6B,C) represent that conjugation across the benzothiazole and quinoline groups does not exist for the S1 minimum, and the excitation from HOMO to LUMO corresponds to charge transfer from the benzothiazole to the quinoline group. The absence of the conjugation between these two groups can be revealed by a bond order analysis. The Mayer’s bond order of the C−C double bond (labeled with ψ in Figure 6A, left) for the twisted structure is only 0.99; hence, this bond is rather characterized as a single bond, as expected. At the S0 minimum structure, the Mayer’s bond order of the bond is reduced from 1.30 to 1.17 upon excitation. It implies that the excitation breaks the conjugation between the benzothiazole and quinoline groups, and the C−C bond labeled with ψ loses double bond character as a result. Consequently, the rotation around ψ would be triggered to alleviate steric repulsion between the two groups. We then optimized the DNA−TOTO complex at the ground state, and its optimized structure is shown in Figure 5A. The RMS deviation of the optimized geometry with respect to the experimental geometry was 1.021 Å. However, we will not discuss this deviation, because the main purpose of this study is the development of (TD-)DFTB/PCM rather than the parametrization. The absorptions calculated for the optimized complex are 480.1 (44th excitation) and 476.6 (48th excitation) nm. The experimental absorption maximum of the complex is observed at 513 nm, so there is a red-shift of 32 nm upon the binding to DNA. The red-shift predicted by TD-DFTB3/PCM is less than 10 nm; hence, the predicted shift is rather small compared to that in the experiment. The oscillator strengths of these excitations for the complex are 0.121 and 0.486, respectively, while those for the free TOTO are 0.627 (477.2 nm) and 1.129 (471.3 nm). These values suggest that complexation may decrease the absorbance, and the observation is qualitatively consistent with the observed experimental absorbance decrease.97 The excitation vectors indicate that the excitations in the complex mainly consist of local excitation in the TOTO dye. However, there are some contributions from environmental DNA base pairs, and this will be discussed later. The excited-state geometry optimizations required an immense effort, since these states easily alter during optimization. Nevertheless, we have managed to optimize the complex on the two excited-state PESs relevant to local TOTO excitations discussed above. Note that we did not consider any energy transfer through DNA or state crossings such as conical intersections. Instead, we simply followed a state with the desired excitation vectors during optimization. We were able to confirm that the excitation character is reasonable by analyzing the MOs and excitation vectors at the initial and final geometries. The calculated emission lies at 539.0 and 614.6 nm, which correspond to the 34th and 16th excitations, respectively, at geometry convergence. The dihedral angle ψ is −3.5° and −29.9°, indicating that the environmental DNA suppresses rotation around ψ. The orbitals relevant to these excitations are shown in Figure 7. For the emission at 539.0 nm, three orbitals exist that contribute predominantly to the excitation. The HOMO−6 → LUMO excitation is the most dominant contribution, with around 50% mixing, and corresponds to local TOTO excitation. The HOMO−22 → LUMO and HOMO−18 → LUMO excitations contribute 24.7% and 14.7% to the mixing, respectively; these occupied

Figure 7. Orbitals [±0.04 (e/a0)1/2 isosurface] relevant to emission. Values in parentheses represent the eigenvalues of vectors in Hartree.

orbitals are on neighboring base pairs. In contrast, the emission at 614.6 nm is 95.2% comprised of local TOTO excitation. Thus, it seems that the emission maximum has a strong dependence on the mixing of orbitals relevant to excitation. However, it is unfortunately difficult to identify which emission corresponds to the experimental emission at 532 nm97 using the current methodology. Nevertheless, considering that the experimental Stokes shift is 19 nm, it is reasonable to deduce that the emission at 539.0 nm corresponds to the experimentally observed emission at 532 nm. The mixing of orbitals from environmental DNA for the emission at 539.0 nm in Figure 7 indicates that the distances between TOTO and base pairs are lesser for the emission at 539.0 nm than those at 614.6 nm. In this case, the result implies that excitation from environmental DNA base pairs has to be included in linearresponse TD calculations to study such a mechanism, which depends on the distances between the dye and the DNA base pairs. We must also comment on the mechanism for the enhancement of emission upon complexation. When TOTO is free, nonradiative decay is dominant as a result of the rotation around the double bond discussed above. Although charge transfer deexcitation or state-crossing may occur, we cannot give a detailed description of the decay mechanism. When complexed, the TOTO molecule can no longer move freely, so emission occurs efficiently. Given that the CAM-B3LYP calculations also found another minimum on the S1 PES with ψ = −40.0°, the minimum may also play a role in the real emission mechanism.

■

CONCLUSIONS We derived an analytical first-order geometrical derivative of ground and excited state energies within the framework of TDDFTB/PCM and implemented it in GAMESS-US. The analytical gradient of the excited state energy required solving Z-vector equations, similar to that for the case without PCM. J


Article


We expect that the method developed in this study will be available as part of GAMESS in the near future. In future work, (TD-)DFTB/PCM will be combined with fragment molecular orbital-based DFTB.54,55 Long-range-corrected DFTB52,53 is also promising. LC-BLYP improved the prediction of absorption and fluorescence energies over BLYP, so longrange-corrected DFTB is also expected to improve the results, as very recently shown in ref 53.

The implemented gradient was found to be fairly accurate in comparison with numerical gradients, and in contrast to the previous implementation,33 our new DFTB339 calculations can be combined with PCM in the developed code. With the latter new tool, we examined the performance of (TD-)DFTB/PCM, focusing on the quality of the calculated PES. At the ground state, a set of halogen-exchange SN2 reactions was tested with DFTB2/PCM and DFTB3-D3(BJ)/ PCM. For isolated CH3X (X = F, Cl, Br, and I), DFTB3D3(BJ) reproduced the C−X bond distances obtained with B3LYP/aug-cc-pVTZ very well, although this result is not that emphatic because the DFTB3 parameters for the halogen atoms were parametrized43 to reproduce the B3LYP/def2TZVPP geometries. For TSs, it also gave a prediction consistent with B3LYP/aug-cc-pVTZ for Cl and I, but did give somewhat inconsistent results for F and Br. In particular, the latter gave a large deviation of 0.096 Å for the C−Br distance and more than 10 kcal/mol for ΔG compared with the B3LYP/aug-cc-pVTZ result. The analytical excitation energy gradient makes the exploration of excited-state PESs efficient. With 3HF, focusing on the ESIPT reaction, TSs for the ground and the first singlet−singlet excited states were obtained. (TD-)DFTB3 reproduced the relative energies obtained with B3LYP/aug-ccpVDZ//B3LYP/6-31G(d,p) within a few kcal/mol. However, the prediction of the experimental absorption and fluorescence energies with TD-DFTB was rather poor. Additional DFT calculations showed that PBE and BLYP functionals also underestimated these energies significantly, and DFTB3 reproduced the PBE and BLYP values within around 0.1 eV. As the electronic part of the DFTB parameters were obtained with the PBE functional, this result is in a sense satisfactory. Using a long-range corrected functional, LC-BLYP, the experimental results were reproduced very well, which indicates that long-range corrections are important for DFTB, as well as DFT. In terms of computational efficiency, TD-DFTB/PCM (single point) calculations for systems consisting of 1000 atoms are manageable. For the smallest system in the example (C100H102; finished in 20 s), the additional cost due to PCM was 1.42 times greater than the TD-DFTB calculation in a vacuum. However, for the biggest system (C500H502; finished in 40 min with one CPU core), the additional cost was only 10%. The calculation for 1IO5 protein shows that the convergence of the ASC equation (eq 8) becomes hard, compared to that for the nonpolar molecule. Even with such an additional cost for TD-DFTB/PCM, it is significantly faster than ordinary TD-DFT/PCM. For instance, a vibrational frequency calculation for 3HF in ethanol with the analytical gradient at TD-B3LYP/6-31G(d,p) took around 28 h with six CPU cores (Xeon E5-1650 v3), but the corresponding calculation with TD-DFTB3/PCM finished in just 5 min. This speed-up, by a factor of 340, with the moderate accuracy of DFTB provides a useful tool for predicting the properties of large systems to which conventional DFT calculations are hard to apply. As a demonstrative example, we performed (TD)DFTB3/PCM calculations for a DNA−TOTO complex (PDB ID 108D) containing 619 atoms. The experimental enhancement of the fluorescence intensity in TOTO is attributed to the suppression of the rotation around the double bond between the quinoline and the carbon between the benzothiazole and quinoline groups.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.5b10732. Accuracy of implemented analytical gradient (Table S1), the effect of dispersion corrections for DFTB2 (Table S2), a comparison of RC−X and ΔG for CH3X (X = F, Cl, Br, and I) in vacuum and in water (Table S3), Oa−H distances of stationary points for 3HF (Table S4), Cartesian coordinates of optimized 3HF, and frontier orbitals of free TOTO (Figure S1) (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +81 (0) 75 711-7894. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Y.N. acknowledges the Fukui Institute for the Fundamental Chemistry Fellowship. This work is supported by JSPS KAKENHI Grant Number 15H06316 (Grant-in-Aid for Research Activity Start-up).

■

REFERENCES

(1) Lin, H.; Truhlar, D. G. QM/MM: What Have We Learned, Where Are We, and Where Do We Go from Here? Theor. Chem. Acc. 2007, 117, 185−199. (2) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical Treatment of Solvation for Molecular Mechanics and Dynamics. J. Am. Chem. Soc. 1990, 112, 6127−6129. (3) Qiu, D.; Shenkin, P. S.; Hollinger, F. P.; Still, W. C. The GB/SA Model for Solvation. A Fast Analytical Method for the Calculation of Approximate Born Radii. J. Phys. Chem. A 1997, 101, 3005−3014. (4) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3094. (5) Klamt, A.; Schüürmann, G. COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and Its Gradient. J. Chem. Soc., Perkin Trans. 2 1993, 799−805. (6) Ten-no, S.; Hirata, F.; Kato, S. A Hybrid Approach for the Solvent Effect on the Electronic Structure of a Solute Based on the RISM and Hartree-Fock Equations. Chem. Phys. Lett. 1993, 214, 391− 396. (7) Sato, H.; Hirata, F.; Kato, S. Analytical Energy Gradient for the Reference Interaction Site Model Multiconfigurational Self-Consistent-Field Method: Application to 1,2-Difluoroethylene in Aqueous Solution. J. Chem. Phys. 1996, 105, 1546−1551. (8) Miertuš, S.; Scrocco, E.; Tomasi, J. Electrostatic Interaction of a Solute with a Continuum. A Direct Utilizaion of Ab Initio Molecular Potentials for the Prevision of Solvent Effects. Chem. Phys. 1981, 55, 117−129.

K


Article

The Journal of Physical Chemistry A (9) Miertuš, S.; Tomasi, J. Approximate Evaluations of the Electrostatic Free Energy and Internal Energy Changes in Solution Processes. Chem. Phys. 1982, 65, 239−245. (10) Cossi, M.; Mennucci, B.; Cammi, R. Analytical First Derivatives of Molecular Surfaces with Respect to Nuclear Coordinates. J. Comput. Chem. 1996, 17, 57−73. (11) Barone, V.; Cossi, M. Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model. J. Phys. Chem. A 1998, 102, 1995−2001. (12) Amovilli, C.; Mennucci, B.; Floris, F. M. MCSCF Study of the SN2 Menshutkin Reaction in Aqueous Solution within the Polarizable Continuum Model. J. Phys. Chem. B 1998, 102, 3023−3028. (13) Mennucci, B.; Cammi, R.; Tomasi, J. Excited States and Solvatochromic Shifts within a Nonequilibrium Solvation Approach: A New Formulation of the Integral Equation Formalism Method at the Self-Consistent Field, Configuration Interaction, and Multiconfiguration Self-Consistent Field Level. J. Chem. Phys. 1998, 109, 2798−2807. (14) Cammi, R.; Mennucci, B.; Tomasi, J. Second-Order MøllerPlesset Analytical Derivatives for the Polarizable Continuum Model Using the Relaxed Density Approach. J. Phys. Chem. A 1999, 103, 9100−9108. (15) Cammi, R. Quantum Cluster Theory for the Polarizable Continuum Model. I. The CCSD Level with Analytical First and Second Derivatives. J. Chem. Phys. 2009, 131, 164104. (16) Cammi, R.; Mennucci, B. Linear Response Theory for the Polarizable Continuum Model. J. Chem. Phys. 1999, 110, 9877−9886. (17) Cossi, M.; Barone, V. Time-Dependent Density Functional Theory for Molecules in Liquid Solutions. J. Chem. Phys. 2001, 115, 4708−4717. (18) Scalmani, G.; Frisch, M. J.; Mennucci, B.; Tomasi, J.; Cammi, R.; Barone, V. Geometries and Properties of Excited States in the Gas Phase and in Solution: Theory and Application of a Time-Dependent Density Functional Theory Polarizable Continuum Model. J. Chem. Phys. 2006, 124, 094107. (19) Liu, J.; Liang, W. Analytical Second Derivatives of Excited-State Energy within the Time-Dependent Density Functional Theory Coupled with a Conductor-Like Polarizable Continuum Model. J. Chem. Phys. 2013, 138, 024101. (20) Pierotti, R. A. A Scaled Particle Theory of Aqueous and Nonaqueous Solutions. Chem. Rev. 1976, 76, 717−726. (21) Floris, F.; Tomasi, J. Evaluation of the Dispersion Contribution to the Solvation Energy. A Simple Computational Model in the Continuum Approximation. J. Comput. Chem. 1989, 10, 616−627. (22) Floris, F. M.; Tomasi, J.; Ahuir, J. L. P. Dispersion and Repulsion Contributions to the Solvation Energy: Refinements to a Simple Computational Model in the Continuum Approximation. J. Comput. Chem. 1991, 12, 784−791. (23) Yokogawa, D.; Sato, H.; Sakaki, S. New Generation of the Reference Interaction Site Model Self-Consistent Field Method: Introduction of Spatial Electron Density Distribution to the Solvation Theory. J. Chem. Phys. 2007, 126, 244504. (24) Wada, T.; Nakano, H.; Sato, H. Solvatochromic Shift of Brooker’s Merocyanine: Hartree-Fock Exchange in Time Dependent Density Functional Calculation and Hydrogen Bonding Effect. J. Chem. Theory Comput. 2014, 10, 4535−4547. (25) Koskinen, P.; Mäkinen, V. Density-Functional Tight-Binding for Beginners. Comput. Mater. Sci. 2009, 47, 237−253. (26) Seifert, G.; Joswig, J.-O. Density-Functional Tight Binding - An Approximate Density-Functional Theory Method. WIREs Comput. Mol. Sci. 2012, 2, 456−465. (27) Gaus, M.; Cui, Q.; Elstner, M. Density Functional Tight Binding: Application to Organic and Biological Molecules. WIREs Comput. Mol. Sci. 2014, 4, 49−61. (28) Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R. Construction of Tight-Binding-Like Potentials on the Basis of Density-Functional Theory: Application to Carbon. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 51, 12947−12957. (29) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge Density-

Functional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260−7268. (30) Lu, Z.; Liu, H.; Elstner, M.; Yang, W. In Reviews of Modern Quantum Chemistry; Sen, K. D., Ed.; World Scientific: Singapore, 2002; Chapter 53, pp 1606−1614. (31) Xie, L.; Liu, H. The Treatment of Solvation by a Generalized Born Model and a Self-Consistent Charge-Density Functional TheoryBased Tight-Binding Method. J. Comput. Chem. 2002, 23, 1404−1415. (32) Hou, G.; Zhu, X.; Cui, Q. An Implicit Solvent Model for SCCDFTB with Charge-Dependent Radii. J. Chem. Theory Comput. 2010, 6, 2303−2314. (33) Barone, V.; Carnimeo, I.; Scalmani, G. Computational Spectroscopy of Large Systems in Solution: The DFTB/PCM and TD-DFTB/PCM Approach. J. Chem. Theory Comput. 2013, 9, 2052− 2071. (34) Caricato, M.; Mennucci, B.; Tomasi, J. Solvent Effects on the Electronic Spectra: An Extension of the Polarizable Continuum Model to the ZINDO Method. J. Phys. Chem. A 2004, 108, 6248−6256. (35) Steinmann, C.; Blædel, K. L.; Christensen, A. S.; Jensen, J. H. Interface of the Polarizable Continuum Model of Solvation with SemiEmpirical Methods in the GAMESS Program. PLoS One 2013, 8, e67725. (36) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. J.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (37) Elstner, M. The SCC-DFTB Method and Its Application to Biological Systems. Theor. Chem. Acc. 2006, 116, 316−325. (38) Yang, Y.; Yu, H.; York, D.; Cui, Q.; Elstner, M. Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method: Third-Order Expansion of the Density Functional Theory Total Energy and Introduction of a Modified Effective Coulomb Interaction. J. Phys. Chem. A 2007, 111, 10861−10873. (39) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the SelfConsistent-Charge Density-Functional Tight-Binding Method (SCCDFTB). J. Chem. Theory Comput. 2011, 7, 931−948. (40) Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. J. Chem. Theory Comput. 2013, 9, 338−354. (41) Gaus, M.; Lu, X.; Elstner, M.; Cui, Q. Parameterization of DFTB3/3OB for Sulfur and Phosphorus for Chemical and Biological Applications. J. Chem. Theory Comput. 2014, 10, 1518−1537. (42) Lu, X.; Gaus, M.; Elstner, M.; Cui, Q. Parametrization of DFTB3/3OB for Magnesium and Zinc for Chemical and Biological Applications. J. Phys. Chem. B 2015, 119, 1062−1082. (43) Kubillus, M.; Kubař, T.; Gaus, M.; Ř ezác,̌ J.; Elstner, M. Parameterization of the DFTB3 Method for Br, Ca, Cl, F, I, K, and Na in Organic and Biological Systems. J. Chem. Theory Comput. 2015, 11, 332−342. (44) Niehaus, T. A.; Suhai, S.; Della Sala, F.; Lugli, P.; Elstner, M.; Seifert, G.; Frauenheim, T. Tight-Binding Approach to TimeDependent Density-Functional Response Theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 63, 085108. (45) Nishimoto, Y. Time-Dependent Density-Functional TightBinding Method with the Third-Order Expansion of Electron Density. J. Chem. Phys. 2015, 134, 094108. (46) Fukui, K. The Path of Chemical Reactions - The IRC Approach. Acc. Chem. Res. 1981, 14, 363−368. (47) Yamaguchi, Y.; Schaefer, H. F., III; Osamura, Y.; Goddard, J. A New Dimension to Quantum Chemistry: Analytical Derivative Methods in Ab Initio Molecular Electronic Structure Theory; Oxford University Press: New York, 1994. (48) Wang, Y.; Li, H. Excited State Geometry of Photoactive Yellow Protein Chromophore: A Combined Conductorlike Polarizable Continuum Model and Time-Dependent Density Functional Study. J. Chem. Phys. 2010, 133, 034108. (49) Heringer, D.; Niehaus, T. A.; Wanko, M.; Frauenheim, T. Analytical Excited State Forces for the Time-Dependent DensityL


Article

The Journal of Physical Chemistry A Functional Tight-Binding Method. J. Comput. Chem. 2007, 28, 2589− 2601. (50) Furche, F.; Ahlrichs, R. Adiabatic Time-Dependent Density Functional Methods for Excited State Properties. J. Chem. Phys. 2002, 117, 7433−7447. (51) Handy, N. C.; Schaefer, H. F., III On the Evaluation of Analytic Energy Derivatives for Correlated Wave Functions. J. Chem. Phys. 1984, 81, 5031−5033. (52) Niehaus, T. A.; Della Sala, F. Range Separated Functionals in the Density Functional Based Tight-Binding Method: Formalism. Phys. Status Solidi B 2012, 249, 237−244. (53) Humeniuk, A.; Mitrić, R. Long-Range Correction for TightBinding TD-DFT. J. Chem. Phys. 2015, 143, 134120. (54) Nishimoto, Y.; Fedorov, D. G.; Irle, S. Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method. J. Chem. Theory Comput. 2014, 10, 4801−4812. (55) Nishimoto, Y.; Fedorov, D. G.; Irle, S. Third-Order DensityFunctional Tight-Binding Combined with the Fragment Molecular Orbital Method. Chem. Phys. Lett. 2015, 636, 90−96. (56) Li, H.; Jensen, J. H. Improving the Efficiency and Convergence of Geometry Optimization with the Polarizable Continuum Model: New Energy Gradients and Molecular Surface Tessellation. J. Comput. Chem. 2004, 25, 1449−1462. (57) Su, P.; Li, H. Continuous and Smooth Potential Energy Surface for Conductorlike Screening Solvation Model Using Fixed Points with Variable Areas. J. Chem. Phys. 2009, 130, 074109. (58) York, D. M.; Karplus, M. A Smooth Solvation Potential Based on the Conductor-Like Screening Model. J. Phys. Chem. A 1999, 103, 11060−11079. (59) Zheng, G.; Witek, H. A.; Bobadova-Parvanova, P.; Irle, S.; Musaev, D. G.; Prabhakar, R.; Morokuma, K.; Lundberg, M.; Elstner, M.; Köhler, C.; et al. Parameter Calibration of Transition-Metal Elements for the Spin-Polarized Self-Consistent-Charge DensityFunctional Tight-Binding (DFTB) Method: Sc, Ti, Fe, Co, and Ni. J. Chem. Theory Comput. 2007, 3, 1349−1367. (60) The ChemCraft Web site; http://www.chemcraftprog.com/ index.html (accessed Jan 20, 2016). (61) Niehaus, T.; Elstner, M.; Frauenheim, T.; Suhai, S. Application of an Approximate Density-Functional Method to Sulfur Containing Compounds. J. Mol. Struct.: THEOCHEM 2001, 541, 185−194. (62) Kubař, T.; Bodrog, Z.; Gaus, M.; Köhler, C.; Aradi, B.; Frauenheim, T.; Elstner, M. Parametrization of the SCC-DFTB Method for Halogens. J. Chem. Theory Comput. 2013, 9, 2939−2949. (63) The DFTB Web site; http://www.dftb.org (accessed May 15, 2015). (64) Li, H.; Pomelli, C. S.; Jensen, J. H. Continuum Solvation of Large Molecules Described by QM/MM: A Semi-Iterative Implementation of the PCM/EFP Interface. Theor. Chem. Acc. 2003, 109, 71−84. (65) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (66) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456−1465. (67) Becke, A. D. A New Mixing of Hartree-Fock and Local Densityfunctional Theories. J. Chem. Phys. 1993, 98, 1372−1377. (68) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (69) Perdew, J. P. Density-Functional Approximation for the Correlation Energy of the Inhomogeneous Electron Gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (70) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (71) Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023.

(72) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron Affinities of the First-Row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796−6806. (73) Woon, D. E.; Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. III. The Atoms Aluminum through Argon. J. Chem. Phys. 1993, 98, 1358−1371. (74) Peterson, K. A.; Figgen, D.; Goll, E.; Stoll, H.; Dolg, M. Systematically Convergent Basis Sets with Relativistic Pseudopotentials. II. Small-Core Pseudopotentials and Correlation Consistent Basis Sets for the Post-d Group 16−18 Elements. J. Chem. Phys. 2003, 119, 11113−11123. (75) Peterson, K. A.; Shepler, B. C.; Figgen, D.; Stoll, H. On the Spectroscopic and Thermochemical Properties of ClO, BrO, IO, and Their Anions. J. Phys. Chem. A 2006, 110, 13877−13883. (76) Stevens, W. J.; Basch, H.; Krauss, M. Compact Effective Potentials and Efficient Shared-Exponent Basis Sets for the First- and Second-Row Atoms. J. Chem. Phys. 1984, 81, 6026−6033. (77) Stevens, W. J.; Krauss, M.; Basch, H.; Jasien, P. G. Relativistic Compact Effective Potentials and Efficient, Shared-Exponent Basis Sets for the Third-, Fourth-, and Fifth-Row Atoms. Can. J. Chem. 1992, 70, 612−630. (78) Feller, D. The Role of Databases in Support of Computational Chemistry Calculations. J. Comput. Chem. 1996, 17, 1571−1586. (79) Schuchardt, K. L.; Didier, B. T.; Elsethagen, T.; Sun, L.; Gurumoorthi, V.; Chase, J.; Li, J.; Windus, T. L. Basis Set Exchange: A Community Database for Computational Sciences. J. Chem. Inf. Model. 2007, 47, 1045−1052. (80) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257−2261. (81) Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chem. Acc. 1973, 28, 213−222. (82) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (83) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (84) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A Long-Range-Corrected Time-Dependent Density Functional Theory. J. Chem. Phys. 2004, 120, 8425−8433. (85) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. Hydrogen Bonding and Stacking Interactions of Nucleic Acid Base Pairs: A Density-Functional-Theory Based Treatment. J. Chem. Phys. 2001, 114, 5149−5155. (86) Tucker, S. C.; Truhlar, D. G. Ab Initio Calculations of the Transition-State Geometry and Vibrational Frequencies of the SN2 Reaction of Cl− with CH3Cl. J. Phys. Chem. 1989, 93, 8138−8142. (87) Parthiban, S.; de Oliveira, G.; Martin, J. M. L. Benchmark Ab Initio Energy Profiles for the Gas-Phase SN2 Reactions Y− + CH3X → CH3Y + X− (X,Y = F,Cl,Br). Validation of Hybrid DFT Methods. J. Phys. Chem. A 2001, 105, 895−904. (88) Zhang, J.; Hase, W. L. Electronic Structure Theory Study of the F− + CH3I → FCH3 + I− Potential Energy Surface. J. Phys. Chem. A 2010, 114, 9635−9643. (89) Sengupta, P. K.; Kasha, M. Excited State Proton-Transfer Spectroscopy of 3-Hydroxyflavone and Quercetin. Chem. Phys. Lett. 1979, 68, 382−385. (90) Ameer-Beg, S.; Ormson, S. M.; Brown, R. G.; Matousek, P.; Towrie, M.; Nibbering, E. T. J.; Foggi, P.; Neuwahl, F. V. R. Ultrafast Measurements of Excited State Intramolecular Proton Transfer (ESIPT) in Room Temperature Solutions of 3-Hydroxyflavone and Derivatives. J. Phys. Chem. A 2001, 105, 3709−3718. (91) Gonzalez, C.; Schlegel, H. B. An Improved Algorithm for Reaction Path Following. J. Chem. Phys. 1989, 90, 2154−2161. (92) Caricato, M.; Mennucci, B.; Tomasi, J.; Ingrosso, F.; Cammi, R.; Corni, S.; Scalmani, G. Formation and Relaxation of Excited States in M


Article

The Journal of Physical Chemistry A Solution: A New Time Dependent Polarizable Continuum Model Based on Time Dependent Density Functional Theory. J. Chem. Phys. 2006, 124, 124520. (93) Improta, R.; Barone, V.; Scalmani, G.; Frisch, M. J. A StateSpecific Polarizable Continuum Model Time Dependent Density Functional Theory Method for Excited State Calculations in Solution. J. Chem. Phys. 2006, 125, 054103. (94) Improta, R.; Scalmani, G.; Frisch, M. J.; Barone, V. Toward Effective and Reliable Fluorescence Energies in Solution by a New State Specific Polarizable Continuum Model Time Dependent Density Functional Theory Approach. J. Chem. Phys. 2007, 127, 074504. (95) Niimura, N.; Minezaki, Y.; Nonaka, T.; Castagna, J.-C.; Cipriani, F.; Høghøj, P.; Lehmann, M. S.; Wilkinson, C. Neutron Laue Diffractometry with an Imaging Plate Provides an Effective Data Collection Regime for Neutron Protein Crystallography. Nat. Struct. Biol. 1997, 4, 909−914. (96) Spielmann, H. P.; Wemmer, D. E.; Jacobsen, J. P. Solution Structure of a DNA Complex with the Fluorescent Bis-Intercalator TOTO Determined by NMR Spectroscopy. Biochemistry 1995, 34, 8542−8553. (97) Rye, H. S.; Yue, S.; Wemmer, D. E.; Quesada, M. A.; Haugland, R. P.; Mathies, R. A.; Glazer, A. N. Stable Fluorescent Complexes of Double-Stranded DNA with Bis-Intercalating Asymmetric Cyanine Dyes: Properties and Applications. Nucleic Acids Res. 1992, 20, 2803− 2812. (98) Benson, S. C.; Mathies, R. A.; Glazer, A. N. Heterodimeric DNA-Binding Dyes Designed for Energy Transfer: Stability and Applications of the DNA Complexes. Nucleic Acids Res. 1993, 21, 5720−5726. (99) Benson, S. C.; Singh, P.; Glazer, A. N. Heterodimeric DNAbinding Dyes Designed for Energy Transfer: Synthesis and Spectroscopic Properties. Nucleic Acids Res. 1993, 21, 5727−5735. (100) Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid ExchangeCorrelation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51−57.

N


PCM Applied to Ground and Excited State Potential Energy

Recommend Documents