Computational Approach for Studying Optical Properties of DNA

Sep 1, 2016 - Morten S. Nørby , Jógvan Magnus Haugaard Olsen , Casper Steinmann , and Jacob Kongsted. Journal of Chemical Theory and Computation ...
0 downloads 0 Views 1MB Size
Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Computational approach for studying optical properties of DNA systems in solution Morten Steen Nørby, Casper Steinmann, Jógvan Magnus Haugaard Olsen, Hui Li, and Jacob Kongsted J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 01 Sep 2016 Downloaded from http://pubs.acs.org on September 2, 2016

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

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

Page 1 of 29

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

Journal of Chemical Theory and Computation

Computational approach for studying optical properties of DNA systems in solution Morten Steen Nørby,∗,† Casper Steinmann,‡ J´ogvan Magnus Haugaard Olsen,† Hui Li,¶ and Jacob Kongsted† †Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, DK-5230 Odense M, Denmark ‡Centre for Computational Chemistry, School of Chemistry, University of Bristol, Bristol BS8 1TS, United Kingdom ¶Department of Chemistry and Nebraska Center for Materials and Nanoscience, University of Nebraska-Lincoln, NE 68588-0304, USA E-mail: [email protected] Abstract In this paper we present a study of the methodological aspects regarding calculations of optical properties for DNA systems in solution. Our computational approach will be built upon a fully polarizable QM/MM/Continuum model within a damped linear response theory framework. In this approach the environment is given a highly advanced description in terms of the electrostatic potential through the polarizable embedding model. Furthermore, bulk solvent effects are included in an efficient manner through a conductor-like screening model. With the aim of reducing the computational cost we develop a set of averaged partial charges and distributed isotropic dipole-dipole polarizabilities for DNA suitable for describing the classical region in ground-state and excited-state calculations. Calculations of the UV-spectrum of the 2-aminopurine optical probe embedded in a DNA double helical structure are presented. We show that

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

inclusion of polarizabilities in the embedding potential stemming from the DNA double helix is of crucial importance, while the water cluster surrounding the DNA system is well represented using a continuum approach.

1

Introduction

Most known for its capabilities of storing genetic material the DNA molecule holds a central place in all life. In addition, DNA also plays an important role in molecular recognition where it may interact with different protein complexes 1 and smaller organic or organometallic complexes. 2,3 Common to the understanding of such interactions is a detailed knowledge of the electronic structure of DNA and the interacting complexes where especially optical properties may aid in their characterization. Hence the interplay between theory and experiment is of great importance, where on one hand experiments can be used to study real-size systems, while on the other hand theory can be used to gain insight at the atomistic level of model systems mimicking the real-size systems used in experiment. It is thereby of interest to define a suitable computational protocol for describing optical properties of chromophores within biological systems. From a theoretical point of view, modeling of optical properties can be achieved within the framework of response theory. 4 The determining factor for the size of the molecular system that is possible to study using either standardor damped response theory 5–8 is the underlying electronic-structure theory from which the response equations rely upon. This means that most biological systems, such as DNA in its natural environment, will be out of reach for a full treatment based on first principles calculations. Practically one option would be to consider a small model of the full molecular system on which calculations would be feasible. However, excited-state properties of DNA have proven to be environment-sensitive. 9–13 Thus a better solution to this issue is to focus on a smaller part of the system, which can be treated using quantum mechanics (QM), and embed it in the remainder of the system which can be represented by a cruder model such as molecular mechanics (MM). Since most experiments refer to dilute aqueous solutions it 2

ACS Paragon Plus Environment

Page 2 of 29

Page 3 of 29

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

Journal of Chemical Theory and Computation

is often necessary to include bulk solvation in the theoretical description. This has typically been included in two different ways, either explicitly, using an atomistic water model, or implicitly, through a continuum model. 14 Since the introduction of polarizable force fields in QM/MM by Warshel and Levitt 15 it has been shown many times that including polarization effects greatly increases the accuracy of embedding methods. Thus we here use the polarizable embedding (PE) model, taking environment effects into account through a classical polarizable potential. 16 The PE model normally requires an a priori determination of the embedding potential for each molecular structure taken into account in the modeling of the optical properties. Since this is typically a relatively large number of structures, e.g. a series of snapshots taken from a molecular dynamics (MD) simulation, it would be beneficial to construct a library consisting of averaged embedding parameters. This would circumvent the necessity to calculate new parameters for each individual molecular structure and thus allow for substantial savings of computational resources. However, it can potentially have an impact on the accuracy. Thus in the present paper we present and validate averaged embedding potentials for DNA to be used in PE-QM calculations of response properties. Furthermore, we have coupled the PE model with the FixSol solvation model 17 to take bulk effects stemming from the solvent (typically water) into account. The FixSol model belongs to the family of conductor-like screening models (COSMO) 18–22 where the medium is considered as a conductor with the dielectric constant ǫ = ∞ with a subsequent scaling of the induced charge density to recover the effect of having a finite ǫ. The FixSol model differs from other COSMO models by introducing a damping in the interactions between closelying surface charges, thereby avoiding near singularities in the resulting equations. Hence our model will be termed PE-QM/FixSol to distinguish it from other QM/MM/Continuum models. 23–31 With the aim of describing optical properties of large molecular systems we have implemented the PE-QM/FixSol model both within a conventional response theory formalism 5,16 and a damped linear response theory formalism. 7,8,32,33 The latter formalism thus opens for the calculation of optical properties in any part of the electromagnetic spec-

3

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

trum and especially regions with a high density of states are easily addressed. By means of this model we consider system-specific details in connection with the computation of optical properties. This will include importance of dynamical aspects regarding different regions of the system. As a first step, the dynamical effects of the water potential will be considered by addressing several snapshots of the DNA system embedded in explicit water molecules and compare this to a situation where water is represented using the FixSol continuum model. Secondly, the dynamical motion of the DNA helix is addressed by calculating the optical properties for several snapshots and comparing this to the optical properties obtained using a single QM/MM geometry-optimized structure. As a test system we will consider the 2-aminopurine (2AP) embedded in a DNA helix. The 2AP unit is an isomer of the natural nucleic acid base adenine. It can synthetically be incorporated into a DNA string by substitution of an adenine unit where it will form a Watson-Crick base pair with thymine. 2AP is widely used as a probe for studying nucleic acid structure and dynamics 34 due to the somewhat larger fluorescent quantum yield compared to adenine and other nucleic acid bases. 34,35 This makes 2AP incorporated into DNA (see Fig. 1) an ideal candidate as a test system for the implemented PE-QM/FixSol model which is capable of predicting the absorption properties in a non-equilibrium linear response formalism. The paper is structured as follows; In Section 2 we briefly outline the methodological aspects of our approach. In Section 3 we give the computational details behind our calculations, while in Section 4 we discuss the accuracy of the averaged embedding parameters and the computational strategy for obtaining absorption properties of 2AP.

2

Methodology

The PE model consists of an advanced classical description, in terms of distributed multipole moments and electric dipole polarizabilities, of the environment surrounding a region described using quantum mechanics. In this work we augment the standard PE model with

4

ACS Paragon Plus Environment

Page 4 of 29

Page 5 of 29

Journal of Chemical Theory and Computation

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

Journal of Chemical Theory and Computation

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

pole moments and polarizabilities are typically located at the atomic sites constituting the environment. For further details upon the construction of vˆes and vˆind operators we refer to Refs. [ 16,33]. The derivation of the coupling between an atomistic (dipole) polarizable model to a continuum solvation model in a quantum mechanical linear response formalism has been described before and our approach will closely follow previous efforts. 23,24,27 As such, the vˆsol operator in Eq. 2 describes contributions due to the continuum solvation, here included via the FixSol model. 17 Multiscale embedding calculations in general require a number of key decisions and steps by the user. Here we provide a brief overview of the steps required. First the system must be split into a core quantum region and a classical region. Second, an appropriate theoretical level is selected for the quantum region. The PE model can currently be combined with HF, 16 DFT, 16 CC, 36,37 MCSCF, 38 MCSCF-srDFT 39 and SOPPA. 40 Third, parameters for the embedding potential must be obtained. The PE model keeps an atomistic description of the classical region and parameters are derived from separate QM calculations. Alternatively, parameters from standard force fields can be used. For details on the calculation of embedding parameters we refer to Sec. 2.1. To achieve an efficient calculation of excited-state properties, the aim is to treat only the chromophore, in our case the 2AP probe, at the QM level, while keeping an atomistic description of the DNA helix since this represents a highly non-homogeneous environment. Solvation effects from bulk water are efficiently included by a continuum model. To this end we have implemented a QM/MM/Continuum-type embedding model as described above. Even if the coupling of an atomistic polarizable model to a continuum model in a linear response formalism has been described before, 23,24 the difference between this work and previous works lies first of all in the choice of embedding methods, i.e. our implementation is based on the PE model including a very accurate account of the environmental electrostatics. Furthermore, we base our calculation of optical properties on a damped linear response formalism. In this context, the calculation of optical properties is determined without resolving the individual excited states, which is highly beneficial when

6

ACS Paragon Plus Environment

Page 6 of 29

Page 7 of 29

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

Journal of Chemical Theory and Computation

addressing systems and regions of the spectrum where the density of states is high. In the electric-dipole approximation, the linear absorption cross section of a randomly oriented molecular sample is proportional to the imaginary part of the frequency-dependent molecular polarizability and the UV/Vis absorption spectrum can be determined by the calculation of σ(ω) =

ω Im [¯ α(ω)] , ε0 c

(3)

where ε0 is the vacuum permittivity, c is the speed of light, ω is the angular frequency of the radiation, and α ¯ denotes the isotropic polarizability. The molecular polarizability may be determined by linear response theory employing a phenomenological damping term resulting in an expression for the complex molecular polarizability as

ααβ (ω) = −hhˆ µα ; µ ˆβ iiω = µ[1]† E[2] − (ω + iγ)S[2] α

−1

[1]

µβ ,

(4)

where E[2] and S[2] are the electronic Hessian and the generalized overlap matrix, 5 respec[1]†

tively, µ ˆα is the electric dipole operator and µα the electric-dipole property gradient along the molecular axis α. 7,8 The specific contributions to this response function due to the embedding part enters directly the E[2] matrix. 33 Furthermore, by setting γ = 0 we recover the standard response function and spectral properties can thus also be obtained by solving the corresponding generalized eigenvalue equation. 5,16 Within the context of QM response theory we further note that when solving for the surface charges, as required by the continuum (COSMO) formulation, in a response formalism, only the fast solvent response should be considered, i.e., using the optical permittivity of the solvent, ǫ∞ . 41,42 The induced dipoles of the classical region are always considered as a fast response.

7

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

2.1

Construction of Embedding Potentials

The parameters for the classical embedding potential of DNA used in the PE model are derived using the principles of the molecular fractionation with conjugate caps (MFCC) method for distributed properties. 43,44 In MFCC, a large, computationally intractable system is separated into chemically reasonable fragments which in turn are capped with groups of atoms from neighboring fragments. Two neighboring caps are then merged to form a socalled concap fragment. The classical parameters for the entire system are then derived from the calculations on the isolated capped fragments and concaps using the MFCC principle. 44 The computational procedure for constructing the embedding potential of a DNA sequence requires several QM calculations on isolated fragments, which potentially can be a computationally demanding task. However, when considering many snapshots the task of creating an embedding potential for each of these quickly becomes the computational bottleneck. In light of this we seek the possibility of generating a set of averaged parameters, isotropic polarizabilities and charges, that can be used to construct the potential for any DNA sequence in line with recent studies on averaged embedding potentials for solvent molecules. 45 The MFCC procedure has been implemented in an in-house Python script designed to automate the construction of protein embedding potentials. 46 This has later been interfaced with FragIt 47 to allow for the fragmentation of any chemical system. FragIt relies on the SMiles ARbitrary Target Specification (SMARTS) language to identify substructures in a chemical system where covalent bonds are cut. The original SMARTS pattern proposed for DNA systems in FragIt 47 was designed for the fragment molecular orbital 48,49 (FMO) and the effective fragment molecular orbital 50 (EFMO) approaches which rely on other means to handle the interface between fragments than MFCC. Initial tests found that this original pattern produced fragments and concaps that were not suitable for the calculations needed to derive the embedding potential parameters with satisfactory accuracy. In this paper we thus propose, and use, a new SMARTS pattern for PE potentials created using the MFCC procedure defined as ([$(POCC)][$(OC1COCC1)]). The resulting fragments derived from 8

ACS Paragon Plus Environment

Page 8 of 29

Page 9 of 29

Journal of Chemical Theory and Computation

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

Journal of Chemical Theory and Computation

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

snapshots for the subsequent property calculations were performed using the Amber14 package 52 utilizing the ff14SB force field. The standard ff14SB force field does not include charges for the 2AP probe but these were obtained from Ref. 53. A total of 24 Na+ counter-ions were added to the system, using default Amber distances, to ensure a zero net charge for the whole system. Subsequently, the system (DNA + ions) was solvated, adding a buffer region of 12 ˚ A (6402 molecules) around the DNA helix, using the TIP3P water model in a truncated octahedral box. Simulations were performed with periodic boundary conditions and electrostatic interactions were modeled using the particle mesh Ewald method 54 utilizing a cutoff distance of 10 ˚ A. The positions of the water and ions were minimized while keeping the DNA strand frozen, followed by a minimization of the entire system. The system was thermally equilibrated using a two-step procedure where the system first was heated from 0K to 300K over 20 ps using the Langevin temperature equilibration scheme (collision frequency of 1.0 ps−1 ) and keeping the volume constant, and then the second equilibration step was continued at 300 K for another 100 ps using constant pressure (1 atm). The SHAKE algorithm 55 was used to fix the hydrogen bonds. After the equilibration phase a 10 ns MD simulation, in the NVT ensemble, was run using a time step of 2 fs from which 50 equidistant snapshots were extracted for electron structure calculations.

3.2

QM/MM Geometry Optimization

The QM region was chosen to consist of the 2AP probe along with the two stacking bases and the bases hydrogen-bonded to these. The part of the DNA backbone connecting the bases was also included in the QM region. The QM and MM regions where capped with hydrogen atoms between sites with entries 161-160, 228-229, 572-571 and 639-640 from the PDB file (first numbers are QM atoms, while the last numbers refer to MM atoms). The geometry optimization was performed utilizing M06/6-31G* 56–58 level of theory for the QM region and the OPLS2005 force field 59 for the MM region, using the QSite program. 60

10

ACS Paragon Plus Environment

Page 10 of 29

Page 11 of 29

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

Journal of Chemical Theory and Computation

3.3

Embedding Parameters

Following previous studies by Olsen et al. 61 all embedding parameters where calculated using the aug-cc-pVDZ basis set. However, due to the size and negative charge of the fragments we use the CAM-B3LYP functional 62 over the B3LYP functional, 63,64 in order to avoid selfinteraction errors. 65 For the purpose of constructing the set of averaged parameters we have chosen four different DNA sequences, all of which have been subject to a MD simulation in accordance with Di Meo et al. 66 Each DNA sequence consists of a mixture between two complementary bases with 40 bases in total. In detail the sequences are double stranded DNA homo-oligonucleotide sequences d(AT)20 , d(GC)20 , d(A)20 ·d(T)20 and d(G)20 ·d(C)20 . From each of these simulations a single snapshot was extracted and the embedding parameters was created according to the MFCC procedure as descriped in Sec. 2.1. The embedding parameters from each of these snapshots where then sorted into four different categories determined by the nucleotide and the averaged parameters were constructed by the mean value of all charges and polarizabilities associated with a unique atom site. However, when constructing the averaged parameters the 8 outermost bases were excluded from each of the considered snapshots in order to prevent edge effects when considering finite DNA sequences. The averaged parameter set consisting of fitted charges and isotropic polarizabilities for each atom-type can be found in the supporting information. Fitted charges where obtained by means of the restrained electrostatic potential (RESP) charge fitting scheme 67 using the Merz-Kollman (MK) scheme 68,69 in Gaussian 09 70 to generate the electrostatic potentials with ten molecular surfaces defined by 1.4 to 2.7 times the vdW surface and a density of 17 grid points per square Angstrom. Fitting of the partial charges was performed using the Antechamber module 71 of the AmberTools15 package. 52 Isotropic polarizabilities where computed using the LoProp procedure 72 using the newly developed implementation based on response functions. 73 Furthermore, when considering geometry-specific embedding potentials, i.e. embedding potentials derived from parameters not subjected to any averaging, distributed atom-centered multipole moment and (an)isotropic polarizabilities were calcu11

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

lated in a similar way as the latter isotropic polarizabilities.

3.4

Optical Properties

The near UV absorption spectrum of the 2AP probe was computed using damped linear response theory combined with the newly developed PE-QM/FixSol model implemented in a development version of the Polarizable Embedding library (PElib) 74 which is interfaced to the Dalton program package. 75 All calculations include only the 2AP probe in the quantum region. All zeroth order multipole moments in a radius of 1.3 ˚ A of the quantum region were redistributed to the three closest sites of the respective expansion site while higher order moments in this region where deleted. The FIXPVA2 tesselation scheme as implemented in the QuanPol package 31 was interfaced with PElib. We used the CAM-B3LYP 62 functional and the pcseg-1 basis set by Jensen 76 unless otherwise noted. When including FixSol in the calculations, atomic radii of 1,40, 2.10, 2.00, 1.90, and 2.20 ˚ A were used for H, C, N, O, and P atoms, respectively, to define the molecular cavity. For the ground-state wave-function optimization we use a static dielectric constant of ǫ = 78.39 whereas for the transition property calculations we used ǫopt = 1.776. When evaluating the response function, Eq. 4, we used a common empirical linewidth factor γ equal to 1000 cm−1 (0.00456 a.u.) and the step between successive frequencies in the transition region equal to 0.0025 a.u.

4

Results and discussion

In the following section we will first consider the development of a set of averaged parameters for a DNA polarizable potential. The central element in electrostatic and polarized embedding calculations is the electrostatic potential calculated from the set of multipoles and polarizabilities because this is what enters directly in the embedding operator (Eq. 2). Therefore consequently our analysis of the quality of the embedding parameters, being either geometry-specific or averaged, will be focused on their ability to reproduce the electrostatic

12

ACS Paragon Plus Environment

Page 12 of 29

Page 13 of 29

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

Journal of Chemical Theory and Computation

potential calculated by a reference QM method as described in ref 61. As a measure for the quality of the classical ESP we calculate the root mean-square deviation (RMSD) of the ESPs compared to a reference ESP on a molecular surface defined by twice the vdW atomic radii for all atoms. The RMSDs are reported in terms of energies from the interaction with a unit point charge. Finally, the set of averaged parameters is used for the computation of optical properties for the embedded 2AP probe. We adopt the notation Mx, where x indicates the level at which the multipole expansion will be truncated. Thus M0 indicates only charges, M1 includes charges and dipoles etc. Charges calculated using the RESP fitting scheme are denoted R0. We use a similar notation for the polarizabilities. These will be denoted as Py, where y = 1 for isotropic polarizabilities and y = 2 for anisotropic polarizabilities. Due to the specific construction of the polarizabilities we exclude polarization between atoms that have been part of the same fragment during the MFCC procedure to avoid double counting of intrafragment polarization. Table 1: Accuracy of the ESP of single nucleotides as a function of multipole order (M0: charges, M1: charges and dipoles, M2: charges, dipoles and quadrupoles, R0: RESP charges) relative to the ESP calculated at CAM-B3LYP/aug-cc-pVDZ level of theory. The ESPs were evaluated on a molecular surface defined by twice the vdW atomic radii. All values are RMSDs given in terms of energies from the interaction with a unit point charge (kJ/mol). Emb. Pot. M01 M11 M21 R02 R03 ff14 1 2 3

Adenine 11.9 23.4 3.0 2.1 8.4 18.6

Thymine 8.1 18.3 2.8 1.7 10.0 14.2

Guanine 9.7 22.9 3.0 2.4 10.6 23.7

geometry-specific LoProp parameters geometry-specific RESP charges averaged RESP charges (cap level 4)

13

ACS Paragon Plus Environment

Cytosine 9.0 19.4 3.0 2.1 10.0 10.7

Journal of Chemical Theory and Computation

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

4.1

Accuracy of the averaged charges

We start by benchmarking the performance of the multipoles and (averaged) fitted charges by examining the ESP created by these parameters for each of the individual nucleotides in isolation (structures have been obtained by extracting a single nucleotide from a snapshot of the MD simulation described in Sec. 3.1 and subsequently performing a geometry optimization utilizing the B3LYP 63,64 functional and the 6-31+G* 57,58,77,78 basis set.). These results are presented in Table 1. The presented numbers provide a direct measure of the errors introduced by using either a multipole expansion or fitted charges compared to using a full quantum-mechanical description when evaluating the ESP. Starting with the ESP generated by the geometry-specific LoProp potentials we observe that the M1 potential is consistently worse than the M0 potential. This spurious convergence pattern with respect to multipole order is in agreement with previous studies. 16,61,79 The RMSDs associated with the LoProp M2 potential are in all cases very small (≈ 3 kJ/mol) which is also consistent with the beforementioned studies. We also analyze the performance of the geometry-specific fitted charges and the averaged fitted charges. For all cases the geometry-specific RESP charges perform slightly better than the parameters from the M2 embedding potential with RMSDs ranging from 1.7 to 2.4 kJ/mol. Furthermore the RMSDs for the averaged charges are comparable to the LoProp M0 potential, ranging from 8.4 to 10.6 kJ/mol. Not surprisingly, we find that charges from the Amber ff94 force field (which are also used in many other variants of Amber force fields such as ff99 and more recent variations of this) are outperformed by the other embedding potentials, including the averaged charges (R0), most likely because these are computed at a lower level of theory (HF/6-31G*). However, we note that the decrease in error when going from Amber charges to averaged R0 charges is substantially smaller for the thymine and cytosine bases.

14

ACS Paragon Plus Environment

Page 14 of 29

Page 15 of 29

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

Journal of Chemical Theory and Computation

4.2

Errors from fragmentation and capping

Analyzing systems with more than one DNA base will allow for a closer examination of the errors stemming from the fragmentation, herein the use of different cap levels and the effect introduced by including polarizabilities in the embedding potential. We start by considering three stacking bases (thymine, guanine and cytosine) connected by sugar and phosphate backbone. Fig. 3 shows the ESP generated by the averaged parameters in terms of RMSDs to the full QM ESP as a function of distance to the van der Waals surface. The RMSDs are averaged over 5 structures, where the structures were extracted from the MD simulation described in Sec. 3.1. We include results obtained using different capping levels to generate the averaged parameters. The large RMSDs are connected to the fact that in order to obtain a full QM reference we must consider a small finite DNA sequence whereas the averaged parameters are developed by excluding edge effects as discussed in Sec. 2.1. At short distances to the van der Waals surface we observe a rapid increase in the RMSD. This is attributed to charge penetration errors. 80 When considering averaged parameters computed using cap level 4 compared to using cap level 2 and 3, we observe a decrease in RMSD values (see Fig. 3) of the same magnitude as inclusion of polarizabilities. Table 2 shows results where we study a trimer model of DNA, consisting of three adeninethymine base pairs connected by the sugar and phosphate backbone. First we notice that when neglecting polarizabilities in the LoProp potential M2, the resulting RMSDs are lowered by 3.1 kJ/mol while decreasing the cap level from four to two. Furthermore, when considering embedding potentials that include polarization effects, the convergence with respect to cap level can be expected to converge faster than for the potentials consisting only of distributed multipole moments. By inspection of Table 2 we also see that potentials based on cap level 2 have smaller errors than potentials based on capping level 4, i.e. for M0P1, M2P2 and M2P1 potentials. However, the differences are rather small and, importantly, much smaller than the error from neglecting polarization. Therefore we have not pursued the behavior of embedding potentials using higher cap levels any further. Finally, we test the averaged 15

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

R0P1 Cap lvl. 2 R0P1 Cap lvl. 3 R0P1 Cap lvl. 4

R0 Cap lvl. 2 R0 Cap lvl. 3 R0 Cap lvl. 4

30 25 RMSD / kJ/mol

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

Page 16 of 29

20 15 10 5 0 1.2

1.4

1.6

1.8 2 2.2 2.4 vdw radius factor

2.6

2.8

Figure 3: Distance dependence of the RMSD of the ESP produced by embedding potentials based on averaged parameters relative to the QM ESP of the three stacking nucleotides, thymine, guanine and cytosine. parameters for this system. Here the effect of adding polarizabilities in the potential more than halves the RMSDs, e.g. for capping level 4 the RMSD goes from 23.3 to 11.2 kJ/mol. In summary, we find that it is necessary to use the geometry-specific LoProp multipoles and polarizabilities, i.e. the M2P2 potential, in order to obtain highly accurate potentials. However, the task to recalculate all parameters for each individual snapshot would be very expensive in the next part of this study where 50 snapshots are extracted from an MD simulation and subsequently used in linear response calculations. This computational expense is avoided by using the set of averaged parameters (in the following we use cap level 2), but with a reduction of accuracy as discussed above.

16

ACS Paragon Plus Environment

Page 17 of 29

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

Journal of Chemical Theory and Computation

Table 2: Accuracy of the trimer model system ESP produced by different potentials relative to CAM-B3LYP/aug-cc-pVDZ evaluated on a molecular surface defined by twice the vdW atomic radii. The ESPs were evaluated on a molecular surface defined by twice the vdW atomic radii. All values are RMSDs given in terms of energies from the interaction with a unit point charge (kJ/mol).

Cap lvl. 2 3 4 1 2

4.3

1

M2 22.0 20.8 18.9

Embedding potential M0P11 M2P21 M2P11 R02 R0P12 11.6 5.1 8.3 25.4 11.5 12.2 5.1 9.4 25.8 12.8 12.6 6.5 10.2 23.3 11.2

geometry-specific LoProp parameters averaged parameters

Optical Properties of 2-Aminopyrine inside DNA

Having examined the errors introduced by using a well-defined set of parameters to describe the interaction of the environment with the quantum region, we now turn the attention to calculations of the optical properties of the 2AP probe embedded in a DNA strand and examine different levels of approximations associated with the computation of these properties. As a reference we will use the near UV absorption spectrum from PE-QM calculations where the 2AP probe is in the QM region, the remaining DNA is described by means of the set of averaged parameters (R0P1). The water molecules are represented with the very accurate M2P2 potential. Since we use the same TIP3P water geometry for each water molecule the parameters for this potential will be the same for each water molecule although rotated according to the specific orientation of each water molecule. Given the errors due to the averaged charges compared to the geometry-specific charges presented in Table 1 an alternative to using averaged charges, one could compute the geometry-specific charges for a single snapshot and further use this set of parameters throughout the averaging over snapshots. To capture dynamical effects we average our results over 50 snapshots from the MD simulation described in Section 3.1. A cumulative average of the absorption profile is constructed using the results obtained from the PE-QM calculation including explicit water molecules. The result, shown in supporting information Figure SI-2, demonstrates that the averaged 17

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

spectrum is converged by including 50 snapshots in the averaging. To a first approximation of the reference system, we remove all explicit water molecules and represent these using the FixSol model. By comparing these results and contrasting to the results obtained with the M2P2 water potential we get an estimate of the importance of explicit solvation. Finally, we consider the approximation of using a single QM/MM geometry-optimized structure to calculate the near UV absorption spectra which neglects all nuclear dynamical effects. In all three approximations we keep the Na+ counter-ions to ensure an overall zero charge and utilize average parameters developed using capping level 2. Fig. 4 shows the computed near UV spectra of the 2AP probe for the three cases described above using the damped response formalism. First we note that in all three spectra we observe three absorption bands. Although, when considering the PE-QM/FixSol combination with no explicit waters, the second band is only seen as a small shoulder to the more intense close-lying absorption band. This is in accordance with previous theoretical and experimental studies of 2AP. 81 Direct comparison of peak positions and intensities to experiment is not reasonable since the experimentally obtained spectrum was that of 2AP in solution and not incorporated in a DNA sequence. A close inspection of the response vectors allows for a complete assignment of the electronic spectrum. In accordance with Holm´en et al. 81 we find the first absorption band, at ≈ 4.3 eV, is due to a π → π ∗ transition (S0 → S1 ). The second band, at ≈ 5.6 eV, we find to be a n → π ∗ (S0 → S2 ), while the third absorption band is due to an intense π → π ∗ transition (S0 → S3 ). Comparing the three theoretically obtained spectra in Fig. 4 we consider the full PE-QM calculation with explicit water as our reference. By removing all explicit water and replacing this by a continuum description through FixSol the position and intensity of the S0 → S1 transition is nearly indistinguishable. Furthermore, a 0.2 eV redshift is observed for the S0 → S3 transition compared to the reference. On the other hand the peak intensities are maintained satisfactorily. When investigating the spectrum based on a single QM/MM geometry-optimized structure, we see that both the S0 → S1 transition and the S0 → S3 transition are shifted by 0.3 eV in comparison to the PE-QM spectrum. Furthermore, the

18

ACS Paragon Plus Environment

Page 18 of 29

Page 19 of 29

Journal of Chemical Theory and Computation

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

Journal of Chemical Theory and Computation

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

5

Conclusion

In this paper we presented a computational study of different aspects connected to the modeling of optical properties using the fully polarizable PE-QM/FixSol model in a linear response TDDFT framework utilizing a non-equilibrium solvation scheme. The present implementation of the PE-QM/FixSol method allows for the calculation of optical linear response properties using either standard- or damped response theory. The former allows for a characterization of excitation energies through the poles of the linear response function and transition moments from the residues at the frequency of the excited state can be obtained. The latter opens for a determination of the absorption spectrum in any given region of the energy spectrum and is the ideal theory when considering energy regions with a high density of states. Embedding calculation of optical properties requires a high-quality embedding potential that can be obtained via ab initio calculations on the molecules and fragments that make up the environment. This is a computationally demanding task especially when considering many MD snapshots. Thus, in order to avoid this task we have developed a set of averaged embedding potential parameters suitable for describing DNA double-helix structures. The set includes fitted RESP charges and isotropic dipole-dipole polarizabilities. By benchmarking this set of parameters against the Amber ff94 charges we obtain better agreement with the reference QM calculation. However, a considerable error is still present when using averaged R0 charges over geometry-specific multipoles up to second order (M2). We show that the inclusion of polarizabilities in the force field is of crucial importance giving a much better agreement with the reference calculation in all cases considered. We illustrate our approach by studying the near UV spectrum of 2AP embedded in a DNA specific helix structure. Especially we focus on how to best represent the full system, i.e. chromophore, DNA helix and water solvent. Taking the PE-QM calculation including explicit water molecules as the reference, we have, as an initial step towards a more computationally efficient model, investigated the possibility of removing all water molecules and representing these with a continuum description. The presented pilot calculations show that 20

ACS Paragon Plus Environment

Page 20 of 29

Page 21 of 29

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

Journal of Chemical Theory and Computation

this does not affect the absorption spectrum substantially.

Acknowledgement Computation/simulation for the work described in this paper was supported by the DeIC National HPC Center, SDU. J.K. thanks the Danish Council for Independent Research (the Sapere Aude program), the Villum Foundation, and the Lundbeck Foundation. C.S. thanks the Danish Council for Independent Research (the Sapere Aude program) for financial support (grant ID: 4181-00370). J.M.H.O. acknowledges financial support from the Danish Council for Independent Research through the Sapere Aude research career program. HL acknowledges the University of Nebraska Holland Computing Center for providing resources in the development of the FixSol method.

Supporting Information Available Details of the atom types, averaged partial charges and polarizabilities (Tables S1-S4). Near UV spectrum (Figure S1) of 2AP embedded in a DNA utilizing capping level 4 parameters. Cumulative average of PE-QM calculated spectrum (Figure S2). This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Widom, J. Q. Rev. Biophys. 2001, 34, 269–324. (2) Metcalfe, C.; Thomas, J. A. Chem. Soc. Rev. 2003, 32, 215. (3) Ortmans, I.; Elias, B.; Kelly, J. M.; Moucheron, C.; Kirsch-DeMesmaeker, A. Dalton Trans. 2004, 668. (4) Norman, P. Phys. Chem. Chem. Phys. 2011, 13, 20519. 21

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

(5) Olsen, J.; rgensen, P. J. J. Chem. Phys. 1985, 82, 3235. (6) Helgaker, T.; Coriani, S.; Jørgensen, P.; Kristensen, K.; Olsen, J.; Ruud, K. Chem. Rev. 2012, 112, 543–631. (7) Norman, P.; Bishop, D. M.; rgen Aa. Jensen, H. J.; Oddershede, J. J. Chem. Phys. 2001, 115, 10323. (8) Norman, P.; Bishop, D. M.; rgen Aa. Jensen, H. J.; Oddershede, J. J. Chem. Phys. 2005, 123, 194103. (9) P´erez, A.; Luque, F. J.; Orozco, M. Acc. Chem. Res. 2012, 45, 196–205. (10) Garrec, J.; Patel, C.; Rothlisberger, U.; Dumont, E. J. Am. Chem. Soc. 2012, 134, 2111–2119. (11) Lavery, R.; Zakrzewska, K.; Beveridge, D.; Bishop, T. C.; Case, D. A.; Cheatham, T.; Dixit, S.; Jayaram, B.; Lankas, F.; Laughton, C.; Maddocks, J. H.; Michon, A.; Osman, R.; Orozco, M.; Perez, A.; Singh, T.; Spackova, N.; Sponer, J. Nucleic Acids Res. 2010, 38, 299–313. (12) P´erez, A.; Lankas, F.; Luque, F. J.; Orozco, M. Nucleic Acids Res. 2008, 36, 2379–2394. (13) Etienne, T.; Very, T.; Perp`ete, E. A.; Monari, A.; Assfeld, X. J. Phys. Chem. B 2013, 117, 4973–4980. (14) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999–3094. (15) Warshel, A.; Levitt, M. J. Mol. Bio. 1976, 103, 227–249. (16) Olsen, J. M.; Aidas, K.; Kongsted, J. J. Chem. Theory Comput. 2010, 6, 3721–3734. (17) Thellamurege, N. M.; Li, H. J. Chem. Phys. 2012, 137, 246101. (18) Klamt, A.; Sch¨ uu ¨rmann, G. J. Chem. Soc., Perkin Trans. 2 1993, 799. 22

ACS Paragon Plus Environment

Page 22 of 29

Page 23 of 29

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

Journal of Chemical Theory and Computation

(19) Stefanovich, E. V.; Truong, T. N. Chem. Phys. Lett. 1995, 244, 65–74. (20) Truong, T. N.; Stefanovich, E. V. Chem. Phys. Lett. 1995, 240, 253–260. (21) Barone, V.; Cossi, M. J. Phys. Chem. A 1998, 102, 1995–2001. (22) Klamt, A. J. Phys. Chem. 1995, 99, 2224–2235. (23) Steindal, A. H.; Ruud, K.; Frediani, L.; Aidas, K.; Kongsted, J. J. Phys. Chem. B 2011, 115, 3027–3037. (24) Caprasecca, S.; Curutchet, C.; Mennucci, B. J. Chem. Theory Comput. 2012, 8, 4462– 4473. (25) Cui, Q. J. Chem. Phys. 2002, 117, 4720. (26) Fedorov, D. G.; Kitaura, K.; Li, H.; Jensen, J. H.; Gordon, M. S. J. Comput. Chem. 2006, 27, 976–985. (27) Li, H.; Gordon, M. S. J. Chem. Phys. 2007, 126, 124112. (28) Lipparini, F.; Cappelli, C.; Barone, V. J. Chem. Theory Comput. 2012, 8, 4153–4165. (29) Li, H. J. Chem. Phys. 2009, 131, 184103. (30) Bandyopadhyay, P.; Gordon, M. S.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 2002, 116, 5023. (31) Thellamurege, N. M.; Si, D.; Cui, F.; Zhu, H.; Lai, R.; Li, H. J. Comput. Chem. 2013, 34, 2816–2833. (32) Kauczor, J.; Jørgensen, P.; Norman, P. J. Chem. Theory Comput. 2011, 7, 1610–1630. (33) Pedersen, M. N.; Hedeg˚ ard, E. D.; Olsen, J. M. H.; Kauczor, J.; Norman, P.; Kongsted, J. J. Chem. Theory Comput. 2014, 10, 1164–1171.

23

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

(34) Jones, A. C.; Neely, R. K. Q. Rev. Biophys. 2015, 48, 244–279. (35) Ward, D. C.; Reich, E.; Stryer, L. J. Bio. Chem. 1969, 244, 1228–1237. (36) Sneskov, K.; Schwabe, T.; Kongsted, J.; Christiansen, O. J. Chem. Phys. 2011, 134, 104108. (37) Schwabe, T.; Sneskov, K.; Olsen, J. M. H.; Kongsted, J.; Christiansen, O.; H¨attig, C. J. Chem. Theory Comput. 2012, 8, 3274–3283. (38) Hedeg˚ ard, E. D.; List, N. H.; Jensen, H. J. A.; Kongsted, J. J. Chem. Phys. 2013, 139, 044101. (39) Hedeg˚ ard, E. D.; Olsen, J. M. H.; Knecht, S.; Kongsted, J.; Jensen, H. J. A. J. Chem. Phys. 2015, 142, 114113. (40) Eriksen, J. J.; Sauer, S. P. A.; Mikkelsen, K. V.; Jensen, H. J. A.; Kongsted, J. J. Comput. Chem. 2012, 33, 2012–2022. (41) Mennucci, B.; Cammi, R.; Tomasi, J. J. Chem. Phys. 1998, 109, 2798. (42) Cammi, R.; Frediani, L.; Mennucci, B.; Ruud, K. J. Chem. Phys. 2003, 119, 5818. (43) Zhang, D. W.; Zhang, J. Z. H. J. Chem. Phys. 2003, 119, 3599. (44) Soderhjelm, P.; Ryde, U. J. Phys. Chem. A 2009, 113, 617–627. (45) Beerepoot, M. T.; Steindal, A. H.; List, N. H.; Kongsted, J.; Olsen, J. M. H. J. Chem. Theory Comput. 2016, (46) Olsen;, J. M. H. Development of Quantum Chemical Methods towards Rationalization and Optimal Design of Photoactive Proteins. 2013; http://dx.doi.org/10.6084/m9. figshare.156852. (47) Steinmann, C.; Ibsen, M. W.; Hansen, A. S.; Jensen, J. H. PLoS ONE 2012, 7, e44480. 24

ACS Paragon Plus Environment

Page 24 of 29

Page 25 of 29

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

Journal of Chemical Theory and Computation

(48) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Chem. Phys. Lett. 1999, 313, 701–706. (49) Fedorov, D. G.; Kitaura, K. J. Phys. Chem. A 2007, 111, 6904–6914. (50) Steinmann, C.; Fedorov, D. G.; Jensen, J. H. J. Phys. Chem. A 2010, 114, 8705–8712. (51) Dallmann, A.; Dehmel, L.; Peters, T.; M¨ ugge, C.; Griesinger, C.; Tuma, J.; Ernsting, N. P. Ang. Chem. Int. Ed. 2010, 49, 5989–5992. (52) Case, D.; Berryman, J.; Betz, R.; Cerutti, D.; Cheatham, T.; III,; Darden, T.; Duke, R.; Giese, T.; Gohlke, H.; Goetz, A.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T.; LeGrand, S.; Li, P.; Luchko, T.; Luo, R.; Madej, B.; Merz, K.; Monard, G.; Needham, P.; Nguyen, H.; Nguyen, H.; Omelyan, I.; Onufriev, A.; Roe, D.; Roitberg, A.; Salomon-Ferrer, R.; Simmerling, C.; Smith, W.; Swails, J.; Walker, R.; Wang, J.; Wolf, R.; Wu, X.; York, D.; Kollman, P. AMBER 1015. University of California, San Francisco. (53) Case, D. A. http://personalpages.manchester.ac.uk/staff/Richard.Bryce/ amber/nuc/D2AP_inf.html, accessed date January 15, 2016. (54) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (55) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. J. Comp. Phys. 1977, 23, 327–341. (56) Zhao, Y.; Truhlar, D. G. Theor. Chem. Account 2007, 120, 215–241. (57) Hehre, W. J. J. Chem. Phys. 1972, 56, 2257. (58) Hariharan, P.; Pople, J. Theor. chim. acta 1973, 28, 213–222. (59) Banks, J. L.; Beard, H. S.; Cao, Y.; Cho, A. E.; Damm, W.; Farid, R.; Felts, A. K.; Halgren, T. A.; Mainz, D. T.; Maple, J. R.; Murphy, R.; Philipp, D. M.; Repasky, M. P.;

25

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Zhang, L. Y.; Berne, B. J.; Friesner, R. A.; Gallicchio, E.; Levy, R. M. J. Comput. Chem. 2005, 26, 1752–1780. (60) Qsite version 5.7 (2011). Schr¨odinger, LLC, New York, NY, 2013. (61) Olsen, J. M. H.; List, N. H.; Kristensen, K.; Kongsted, J. J. Chem. Theory Comput. 2015, 11, 1832–1842. (62) Yanai, T.; Tew, D. P.; Handy, N. C. Chem. Phys. Lett. 2004, 393, 51–57. (63) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (64) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623–11627. (65) Jakobsen, S.; Kristensen, K.; Jensen, F. J. Chem. Theory Comput. 2013, 9, 3978–3985. (66) Meo, F. D.; Pedersen, M. N.; Rubio-Magnieto, J.; Surin, M.; Linares, M.; Norman, P. J. Phys. Chem. Lett. 2015, 6, 355–359. (67) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Chem. Phys. 1993, 97, 10269– 10280. (68) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129–145. (69) Besler, B. H.; Merz, K. M.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431–439. (70) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.;

26

ACS Paragon Plus Environment

Page 26 of 29

Page 27 of 29

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

Journal of Chemical Theory and Computation

Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, V.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian09 Revision D.01. Gaussian Inc. Wallingford CT 2009. (71) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. J. Mol. Graph. Model. 2006, 25, 247–260. (72) Gagliardi, L.; Lindh, R.; Karlstr¨om, G. J. Chem. Phys. 2004, 121, 4494. (73) Vahtras, O. LoProp for Dalton. 2014; http://dx.doi.org/10.5281/zenodo.13276. (74) Olsen, J. M. H. PElib: The Polarizable Embedding library (development version). 2014. (75) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekstr¨om, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fern´andez, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; H¨attig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenaes, E.; Høst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jans´ık, B.; Jensen, H. J. A.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjaergaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnaes, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V. V.; Salek, P.; Samson, C. C. M.; de Mer´as, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; Sylvester-

27

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Hvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; ˚ agren, H. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2013, 4, 269–284. (76) Jensen, F. J. Chem. Theory Comput. 2014, 10, 1074–1085. (77) Francl, M. M. J. Chem. Phys. 1982, 77, 3654. (78) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R. J. Comput. Chem. 1983, 4, 294–391. (79) Schwabe, T.; Olsen, J. M. H.; Sneskov, K.; Kongsted, J.; Christiansen, O. J. Chem. Theory Comput. 2011, 7, 2209–2217. (80) Stone, A. J. The Theory of Intermolecular Forces, 2nd ed.; Oxford University Press, UK, 2013. (81) Holm´en, A.; Nord´en, B.; Albinsson, B. J. Am. Chem. Soc. 1997, 119, 3114–3121.

28

ACS Paragon Plus Environment

Page 28 of 29

Page 29 of 29

Journal of Chemical Theory and Computation

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