Using Density Functional Theory To Study Neutral and Ionized

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Using Density Functional Theory to Study Neutral and Ionized Stacked Thymine Dimers #ura Nakarada, Mihajlo Etinski, and Milena Petkovic J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b06493 • Publication Date (Web): 14 Sep 2016 Downloaded from http://pubs.acs.org on September 15, 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.

The Journal of Physical Chemistry A 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 31

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

The Journal of Physical Chemistry

Using Density Functional Theory to Study Neutral and Ionized Stacked Thymine Dimers Ðura Nakarada, Mihajlo Etinski, and Milena Petković∗ Faculty of Physical Chemistry, University of Belgrade Studentski trg 12-16, 11 158 Belgrade, Serbia E-mail: [email protected]

Phone: +381 (0)11 2187 133. Fax: +381 (0)11 2187 133

∗ To

whom correspondence should be addressed

1 ACS Paragon Plus Environment


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

Abstract Stacking interactions in thymine dimers are studied with Density Functional Theory. According to our calculations, six dimers of comparable stability can be prepared at low temperature, but dimerization is impossible at room temperature due to large entropy contribution that accompanies it. Analysis of vibrational anharmonic coupling terms shows that each of the dimers exhibits distinct vibrational dynamics. Properties of electron density in the intermolecular region are used to analyze neutral stacked species and their ionized forms. Bond paths and critical points in the intermolecular region are identified, but a simple relationship between binding energy and total electron density in the intermolecular critical points could not be found due to an uneven electron distribution in the binding region. Reduced density gradient was confirmed to be a useful tool for analysis of weak stacking interactions. Those interactions also affect vertical and adiabatic ionization energies, which are computed to be slightly lower for the dimers compared to the monomer.


Page 2 of 31

Page 3 of 31

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


Introduction π stacking is crucial for functioning of biological systems: it is essential for molecular recognition processes, but also for structure, stability and dynamics of biological macromolecules (DNA, RNA, various proteins), which prompted their extensive analysis in the last two decades. 1–16 Stacking interactions are also important in material science. For example, they are significant for stability of nanostructres formed by discotic liquid crystals which are used in organic electronics, 17 or for host-guest chemistry. 18–21 These large systems are particularly intriguing since their properties are often governed by weak interactions which might be ignored in small systems. 22 Weakness of these non-covalent interactions represents a difficulty both for their experimental and theoretical investigation. 23 A golden standard to theoretically estimate stability of stacked systems is coupled cluster method that employs full treatment of single and double excitations, as well as perturbative treatment of triple excitations, extrapolated to the complete basis set limit, CCSD(T)/CBS. 12,24,25 Although it represents the most reliable method for calculation of interaction energies, it is quite demanding. Significantly less demanding approaches include Møller-Plesset perturbation theory, MP2, 2,3,5–7,13,26,27 as well as Density Functional Theory, DFT. The later approach either relies on functionals with dispersion corrections, 4,17,20,21 or is based on new functionals that account for dispersion effects. 14,28 DFT methods are employed for analysis of large systems which cannot be studied with wavefunction based methods. This was demonstrated by Kim and coworkers who studied adsorption of nucleobases on graphene, as well as electronic and magnetic properties of graphene based devices. 20 DFT is also used for description of large amplitude motion, like construction of potential energy surfaces to study energy changes that arise upon monomer displacements, or investigation of mode couplings. 29 Beside parameters that describe the system in the ground electronic state, excited state properties of π stacked systems are important since they reveal behavior of a system that is subject to perturbations. For example, reliable estimation of ionization energies and electron affinities is a prerequisite to understand mechanisms of radiation induced and charge transfer processes in 3 ACS Paragon Plus Environment


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

living organisms, although appropriate methods must also be searched for in order to study system’s dynamics in electronically excited states. Most of these calculations employ sophisticated wavefunction based approaches. Serrano-Andrés and coworkers 30 demonstrated that high-level quantum chemical methods, like multiconfigurational perturbation methods CASPT2(IPEA) are required even for analysis of nucleobase monomers, which implies necessity of employment of high level methods for analysis of more complex systems. Krylov and coworkers, for example, determined ionization energies of uracil stacked dimers with equation-of-motion coupled-cluster method for ionization potentials with single and double substitutions, EOM-IP-CCSD. However, certain density functionals are able to correctly describe ionized states of stacked species. Krylov and coworkers 31 analyzed a small cluster composed of hydrogen bonded and stacked nucleic bases (stacked adenine-thymine hydrogen bonded dimer) and showed that M06-2X functional can qualitatively describe ionized state of this tetramer, whereas ω B97XD can be used both for qualitative and quantitative determinations. Furthermore, Sun and coworkers 32 recently demonstrated that long-range corrected density functionals 33 can provide reliable results for ionization potentials and electron affinities, comparable to EOM-IP-CCSD values at significantly lower computational cost. It was also shown 34 that their time-dependent variant yields results in good agreement with resolution-of-the-identity coupled cluster (RI-CC2) and EOM-CCSD(T) results. Initially, electrostatic forces were held responsible for π -π interactions, 35,36 until importance of dispersion interactions was recognized. 37 Nowadays, this dilemma is solved with energy component analysis, like symmetry-adapted perturbation theory SAPT, which enables estimation of different energy contributions: attractive (electrostatic, dispersive, inductive) and repulsive (electrostatic and exchange repulsion). Kim and Hobza 24 with their coworkers used this method to study stability of conformers of the benzene dimer and its substituted analogues. Sherrill and coworkers 9 performed a detailed analysis of stability of stacked nucleic acids and found that dispersion interactions are dominant, electrostatic interactions are weaker, whereas induction interactions have the smallest contribution. Sherrill 11 also used this approach to study caton-π interactions and substituent and heteroatom effects on π -π interactions, as well as intercalation of proflavine


Page 4 of 31

Page 5 of 31

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


into a trinucleotide DNA duplex. He demonstrated that dispersion terms are extremely important for stability (and consequently geometry) of these systems, but also emphasized significance of charge–penetration that has stabilizing effect at short distances. Estimation of energy contributions can also be performed with DFT version of the SAPT method, DFT-SAPT, as demonstrated by Riley and Hobza. 12 Kim and coworkers 25 also used DFT variant of SAPT to compare aromatic

π interactions, aliphatic π interactions and non-π stacking interactions. They also concluded that aromatic systems display significant dispersion and electrostatic interactions, and that these large attractive forces between aromatic species cause tight packing in crystal structures. In addition to SAPT, deeper understanding of π stacking is provided by Quantum Theory of Atoms in Molecules. It enables identification of parts of molecules that are subject to these interactions, 17 and also provides with quantitative description of charge transfer within and between specific units. 16 Moreover, empirical expressions for dependence of interaction energy on electron densities and their Laplacians in bond, ring or cage critical points 26,38 for certain stacked species have been established. Reduced density gradient 39 is also a useful tool for analysis of non-covalent interactions, since it enables to identify presence of weak attractive and repulsive interactions which is particularly important to understand processes in biological systems on a molecular level, but also in material science. As one of the building blocks of DNA strands, thymine is an example of a molecule that undergoes π · · · π interactions. For that reason, we have chosen this system as a case study for analysis of stacking interactions. It has already been demonstrated that derivatives of DNA base pairs with substituted benzyl hydrogen atoms generally exhibit stronger stacking interactions compared to the non-substituted analogues, 9,10,13 except in cytosine where methylation (1-methylcytosine) weakens intermolecular interactions. 14 In this contribution we analyze thymine (T) stacked dimers (neutral and ionized forms), as well as influence of methylation on binding energies through investigation of neutral 1-Methyl-thymine (1mT) homodimers. Ability to form thymine dimers at different temperatures is addressed and also possibility of dimers’ dissociation upon temperature increase. The paper is organized as follows: in the next section a brief overview of theoretical



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

Page 6 of 31

methods is presented. The following section contains analysis of optimized structures, computed interaction energies and investigation of dimer’s stability when heated. In the next section ionized species are studied – binding energies, ionization energies, and qualitative discussion of interaction between two monomers through reduced gradient of electron density is presented. The paper is summarized in the final section.

Theoretical methods All quantum chemical calculations were undertaken with the Gaussian09 program package. 40 M06-2X functional 41 was employed in conjunction with the Popple 6-311++G(2df,2pd) basis set. 42–45 M06-2X is a hybrid meta functional designed by Zhao and Thrular which includes parameters that properly describe short range and intermediate range dispersion interactions. It was proven that it yields reliable results for systems with stacking interactions, 46,47 whereas the chosen triple-zeta basis set with diffuse and polarization functions should be flexible enough to enable proper description of weak intermolecular interactions. In order to test performance of different functionlas, thymine dimers were also optimized with ω B97XD, 48 B3LYP, 49,50 CAM-B3LYP 51 and PBE0. 52–54 Calculations with the last three functionals were performed with and without dispersion corrections. All structures were optimized with Tight convergence criteria and an ultrafine integration grid. Reported binding energies capture electrostatic and dispersion interactions corrected for zero point energies, while neglecting thermal contributions: el,BSSE el ZPE ZPE − 2 × Emonomer + Emonomer ∆E = Edimer + Edimer .

(1)

Dimers’ energies include correction for basis set superposition error (BSSE) estimated through counterpoise correction. 40,55 Gibbs free energies and enthalpies were computed by adding corresponding thermal contributions (computed according to the procedure described in the White Papers on the official Gaussian website) to the BSSE corrected electronic energies. Vibrational properties were analyzed by considering normal mode picture. Anharmonic cubic 6 ACS Paragon Plus Environment

Page 7 of 31

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


and quartic force constants were computed from analytic second derivatives by employing finite difference procedure. 56,57 Analysis of electron density was performed with the AIMAll program package. 58

Results and discussion Binding energies Initial structures for thymine dimers (TD) were taken from reference 10. Optimized structures at M06-2X/6-311++G(2df,2pd) level are presented in Figure 1, and their Cartesian coordinates are given in the Supporting Information, Tables S1-S6. As mentioned in the Introduction, the most precise values for binding energies are considered to be those obtained at the CCSD(T)/CBS level. However, it has been shown that certain density functionals are able to provide with reliable qualitative results. Gu et al. 46 showed that M06-2X/DZP++ interaction energies for uracil dimers are in good agreement with CCSD(T)/aug-cc-pVDZ interaction energies and thus used this approach to estimate stability of thymine dimers. Van Mourik 47 also analyzed stability of thymine dimers and concluded that M06-2X appears to be a very promising method to describe stacking interac-

MTD1

MTD1

MTD2

MTD2

MTD3

MTD4

MTD3

MTD4

MTD5

MTD5

MTD6

MTD6

Figure 1: Optimized structures of stacked thymine and 1-methylthymine homodimers. 7 ACS Paragon Plus Environment


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

Page 8 of 31

tions. She also noted that this functional does slightly underestimate binding energy. Since we are not interested in precise determination of interaction energies, we used M06-2X functional to compare relative stability of stacked thymine dimers only to establish whether they do exist under certain conditions. We thus named the dimers according to their binding energies obtained with the M06-2X functional. Although the presented results are not as precise as coupled cluster methods at complete basis set limit would be, they should provide a reliable qualitative picture. In order to compare results obtained with different density functiounals, we computed binding energies with a set of other functionals. Due to the lack of dispersion terms, PBE0, B3LYP and CAM-B3LYP functionals were not able to predict stable stacked dimers and structures converged to hydrogen bonded species. Adding the D3 version of Grimme’s dispersion as implemented in the Gaussian suite of programs 59 enables optimization of stacked species. Stacked dimers were also optimized with ω B97XD functional, which successfully optimized the six structures due to the fact that it includes dispersion terms. Note that all optimizations were performed with Tight optimization criteria, apart from TD1 with the ω B97XD functional – due to convergence problems, in this particular case default optimization criteria were used. Binding energies obtained according to Equation 1 are presented in Table 1. The results obtained with different functionals are comparable. In all cases the binding energies roughly cover the range 30-40 kJ/mol. An exception is TD6 dimer optimized with PBE0-D3 that is more stable than either of its counterparts. This was Table 1: Binding energies in kJ/mol computed with a series of density functionals combined with the 6-311++G(2df,2pd) basis set.

ω B97XD B3LYP-D3 CAM-B3LYP-D3 PBE0-D3

dimer

M06-2X

TD1

-41

-40

-40

-38

-37

TD2

-39

-40

-41

-39

-39

TD3

-38

-36

-37

-35

-34

TD4

-38

-36

-37

-35

-34

TD5

-33

-38

-37

-35

-35

TD6

-30

-31

-33

-33

-53


Page 9 of 31

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


the only example of a failed optimization, since initial stacked structure transformed into a hydrogen bonded one. ω B97XD/M06-2X predicts TD5 to be more/less stable than TD3 and TD4 by 2/5 kJ/mol. B3LYP-D3 and CAM-B3LYP-D3 on the other hand consider all three dimers equally stable, and also TD2 to be 1 kJ/mol below TD1. To summarize, M06-2X, ω B97XD, B3LYP-D3 and CAM-B3LYP-D3 yield comparable interaction energies for the six thymine dimers with slight differences in their stability. PBE0-D3 failed in optimization of one structure.

Dimerization at low and high temperatures Optimized structures of 1-methylthymine dimers (MTD) dimers are displayed in Figure 1 (initial geometries are taken from reference 29, and Cartesian coordinates of the optimized structures are compiled in the Supporting Information, Tables S7-S12). Binding energies of TD and MTD cover ranges -41 to -30 and -50 to -33 kJ/mol, respectively. These results show that methylation does not alter order of stability, and confirm that it strengthens interactions. Accordingly, the most/least stable MTD/TD dimer is more/less stable than either of TD/MTD dimers. The differences between the two sets of results include the following observations: (i) energy differences between dimers TD1-TD4 are small, up to 3 kJ/mol, whereas MTD1 is by 8 kJ/mol more stable than MTD2-MTD5; (ii) TD5 is by 8 kJ/mol less stable than TD1, while MTD5 is equally stable as MTD2-MTD4. Changes in selected thermodynamic properties that result from dimerization at different temperatures are reported in Table 2. Note the decrease in stabilization that follows dimerization at 30K 30 K compared to the value at absolute zero: since ∆G0K dim and ∆Hdim are comparable within the

error limits of the chosen method, this difference is due to change in entropy contributions. At room temperature, T ∆Sdim term overcompensates the enthalpy term, and dimerization ceases to be spontaneous. We have recently analyzed MTD dimers 29 and have demonstrated differences in their relative stability at low and high temperatures. That is, we discussed changes in stability that occur when samples prepared at low temperature are heated. The results presented in Table 2 of the current manuscript deal with dimerization that takes place at 30 K and 298 K, showing that TD 9 ACS Paragon Plus Environment


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

Page 10 of 31

Table 2: Changes of Gibbs free energies ∆Gdim in kJ/mol, enthalpies ∆Hdim in kJ/mol and entropies ∆Sdim in J/(mol·K) that accompany dimerization. All values are computed at M062X/6-311++G(2df,2pd) level at different temperatures (0 K, 30 K and 298 K) and are corrected for BSSE errors. dimer

0K

30 K

298 K

∆Gdim

∆Gdim

∆Hdim

∆Sdim

∆Gdim

∆Hdim

∆Sdim

TD1

-41

-36

-42

-200

11

-39

-168

TD2

-39

-35

-40

-167

9

-37

-154

TD3

-38

-33

-39

-200

13

-36

-164

TD4

-38

-33

-39

-200

13

-36

-164

TD5

-33

-29

-34

-167

16

-31

-158

TD6

-30

-25

-31

-200

17

-28

-151

MTD1

-50

-45

-50

-167

1

-48

-164

MTD2

-42

-37

-42

-167

11

-40

-171

MTD3

-42

-38

-43

-167

10

-41

-171

MTD4

-42

-38

-43

-167

10

-41

-171

MTD5

-42

-37

-43

-200

12

-41

-178

MTD6

-33

-28

-33

-167

15

-30

-151

and MTD species cannot be prepared at room temperature. It must be kept in mind that results presented in Table 2 lack anharmonic contributions: rotationalvibrational coupling, intrinsic anharmonicity and mode couplings. While ro-vibrational coupling can safely be neglected for high frequency vibrational modes, it might affect dynamics of low frequency modes. Additionally, low frequency motion usually takes place on a relatively flat potential energy hyper-surface and significantly deviates from harmonic motion. Concerning high frequency modes, we have already shown that N-H stretching motion in stacked methyl thymine dimers is strongly anharmonic 29 and precise determination of corresponding vibrational frequencies requires construction of high-dimensional models. Anharmonic contributions can be estimated by adding corresponding anharmonic terms to the harmonic Hamiltonian, and redefining the partition functions used for calculation of all thermodynamic quantities. Kim et al. 60 applied this procedure on


Page 11 of 31

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


water dimer and determined change in decrease of entropy that accompanies dimerization of two water molecules when anharmonic corrections are included. They have shown that anharmonicity causes smaller entropy decrease, i.e. it reduces entropic effect. In large systems with large moments of inertia, like thymine dimers, ro-vibrational coupling might most likely be neglected. If that is the case, anharmonic vibrational frequencies might be used to redefine characteristic vibrational temperatures used to estimate ∆Sdim .

Dimers’ fate upon heating Although it is demonstrated that formation of stacked TD and MTD dimers does not take place at 300 K, the question concerning dimers’ fate upon temperature increase remains open, since it depends on the systems’ response to heating. Namely, all dimers are characterized with a low frequency mode νdim that modulates the distance between the two monomers, which might be held responsible for monomers’ separation provided its overtones are significantly populated. Normal mode displacement vectors for νdim in TD1 are displayed in Figure 2: it is obvious that this type of motion increases the distance between the two units. We have identified νdim modes for all TD dimers, and their harmonic frequencies are presented in Table 3. All of them are close to 100 cm−1 . Increase in temperature would cause population of its overtones, which would ultimately cause dimer dissociation. This simplified picture, nevertheless, lacks description of possible energy exchange with other vibrational degrees of freedom. In order to find out whether intramolecular

Figure 2: Dimer TD1: low frequency mode νdim that modifies the distance between two T molecules (for sake of better visibility, normal mode displacements are scaled).



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

Table 3: TD dimers: νdim harmonic frequencies (in cm−1 ) and absolute values of their cubic and quartic anharmonic terms (in cm−1 ) larger than 300 cm−1 . The following labels are used: τMe - rotation of both CH3 groups in TD1 and TD6, and rotation of a single methyl 1 - N-H out-of-plane group in TD2 and TD4; νCH - C-H stretching within a CH3 group; γNH bending. dimer

νdim

TD1

108

anharmonic terms Kνdim ,νdim ,τ 1

= 728

Kνdim ,νdim ,τ 2

= 430

Kνdim ,τMe ,νCH

= 586

Kνdim ,νdim ,τMe ,τMe

= 353

1 Me ,τMe 2 Me ,τMe

TD2

103

TD3

102

/

TD4

102

Kνdim ,νdim ,τ 1

= 403

Kνdim ,νdim ,τ 2

= 465

Kνdim ,νdim ,γ 1

= 345

Kνdim ,νdim ,γ 2

2 NH ,γNH

= 347

Kνdim ,νdim ,τMe ,τMe

= 425

Kνdim ,νdim ,γ 1

= 494

Kνdim ,νdim ,γ 2

= 610

1 Me ,τMe 2 Me ,τMe

TD5

TD6

99

110

1 NH ,γNH

1 NH ,γNH 2 NH ,γNH

vibrational redistribution is pronounced in these dimers, we computed cubic and quartic coupling terms. Those lowest anharmonic terms with absolute values larger than 300 cm−1 are presented in Table 3. Let us first analyze the most stable dimer, TD1. There are no significant cubic force fields that involve νdim , whereas there are two quartic force constants that couple this mode to two vibrational modes that represent methyl rotation. These large force constants indicate involvement of τMe s in νdim ’s dynamics, which might be reflected in energy exchage among them, or in the position of νdim ’s vibrational levels. Provided energy redistribution is pronounced it might either supress or promote dissiciation, which can only be estimated through deeper analysis of system’s dynamics. Concerning other dimers, each of them exhibits distinct behavior. Mode of interest in TD2 is in addition to τMe coupled to a C-H stretching mode, whereas in TD3 it contains no significant anharmonic terms. νdim is in TD4 coupled to two low frequency modes, each of them involves rotation 12 ACS Paragon Plus Environment

Page 12 of 31

Page 13 of 31

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


of a single methyl group. In TD5 it is coupled to N-H out-of-plane bending modes of the N-H· · · O hydrogen bonded fragments, while in TD6 it interacts strongly to τMe in addition to γNH . The reason for diverse couplings in TD dimers lies in existence of different interactions among molecular fragments: provided the two methyl groups are subject to dispersion interactions (TD1, TD2) or are in equivalent positions (TD6), they most likely affect dimer separation; in case N-H· · · O hydrogen bonds are formed, N-H bending motion influences dissociation. From the above discussion, it is obvious that TD3 displays the simplest dynamics concerning νdim , since heating only causes its excitation. 298 K roughly corresponds to 200 cm−1 , which according to the equipartition theorem within harmonic approximation enables population of its second overtone and would most likely not cause dissociation. One must keep in mind that this is an oversimplification since dissociation cannot be described within a normal mode picture, and also harmonic approximation is far too crude. We, nevertheless, used it to test whether energy necessary for the separation of the two units can be gained or lost through coupling with other degrees of freedom. It is obvious from Table 3 that mode interactions cannot be ignored in five dimers. However, based on the values of anharmonic terms that describe coupling between the low frequency mode νdim with other vibrational modes it is not possible to deduce impact of mode interaction on dimer’s dissociation. Namely, energy can be: (i) transferred to νdim causing dissociation, or (ii) taken away from it via modes presented in νdim Table 3 and carried over to other modes according to the Tier model. 61 The above discussion merely implies complexity of interactions between thymine molecules and necessity for more detailed investigations in order to answer the question whether TDs exist at room temperature. Due to the presence of an additional methyl group in MTD, it is expected that their νdim s display even stronger mode couplings and consequently more complicated vibrational dynamics. Since analysis of anharmonic terms indicates complexity of vibrational dynamics but cannot give information whether thermal energy can cause its dissociation, we turned to Born-Oppenheimer molecular dynamics. We took M06-2X/6-311++G(2df,2pd) optimized TD1 geometry and performed calculation at M06-2x/6-31+G(d,p) level. In addition to taking a smaller basis set, com-



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

putational cost was decreased by performing Hessian calculation only at the beginning, and it was subsequently updated through energy gradients. Stepsize was equivalent to 0.5 fs. The temperature was set to 300 K and the propagation lasted 1.5 ps. At the end of the run the structure corresponded to a stable dimer, which indicates that heating the sample with TD1 molecules to 300 K would most likely not cause dissociation.

Topology of electron density A deeper insight into interactions between two monomers in stacked dimers can give Quantum Theory of Atoms in Molecules through analysis of electron density ρ 62,63 in the space between the two monomers (in the following it will be referred as intermolecular region). Electron density redistributes upon dimer formation, which is reflected in appearance of bond paths (BPs, lines that link nearby nuclei which are characterized with maximum electron density) and critical points (CPs, points at which gradient of electron density is equal to zero). Since we are interested in describing non-covalent interactions, we will focus on BPs and CPs (bond critical points, BCPs, ring critical points, RCPs and cage critical points, CCPs) among atoms that do not belong to the same monomer. As indicated by Bader, 64 existence of bond paths confirms presence of attractive interactions. Bond paths and critical points present in TD and MTD species are depicted in Figure 3, and numbers of CPs associated with non-covalent interactions are tabulated in Table 4 (note that the number of bond critical points is identical to the number of bond paths). Generally, bond paths are not straight lines, but are slightly curved towards dimers’ interior. BCPs between heavy atoms (C, N, O) are roughly equally spaced from the two atoms, whereas BCPs that involve H and C/N/O are closer to the hydrogen atom due to the absence of core electrons in the hydrogen atom. Depending on the atom arrangement, number of RCPs is equal or larger than the number of BCPs. Each dimer is characterized by at least one CCP, which is a clear indication of stacking interactions. 26 Due to the presence of an additional methyl group in each monomer, MTDs have more BPs and CPs compared to TD analogues. Although the number of CPs varies from dimer to dimer, larger number of BPs and CPs does not indicate stronger interactions. However, MTDs are 14 ACS Paragon Plus Environment

Page 14 of 31

Page 15 of 31

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


Bond critical points

MTD1

MTD1

MTD2

MTD2

MTD3

MTD4

MTD3

MTD4

MTD5

MTD5

MTD6

MTD6

Ring critical points

MTD1

MTD1

MTD2

MTD2

MTD3

MTD4

MTD3

MTD4

MTD5

MTD6

MTD5

MTD6

MTD5

MTD6

Cage critical points

MTD1

MTD1

MTD2

MTD2

MTD3

MTD4

MTD3

MTD4

MTD5

MTD6

Figure 3: TD and MTD bond critical points, ring critical points and cage critical points. generally characterized with a larger number of bond paths and critical points.



Table 4: Properties of electron density in thymine and 1-methyl-thymine stacked dimers. Labels #BPs and #CP stand for number of bond paths and number of critical points, respectively, that result from dimer formation. Sums of electron density in the bond critical points (∑ ρBCP ), ring critical points (∑ ρRCP ) cage critical points (∑ ρCCP ) and all critical points(∑ ρCP ) between two monomers are given in atomic units. The corresponding energy densities are reported in brackets (also in atomic units).

16

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

Page 16 of 31

dimer

#BCPs

∑ ρBCP (∑ HBCP )

#RCPs

∑ ρRCP (∑ HRCP )

#CCPs

∑ ρCCP

∑ ρCP (∑ HCCP )

TD1

5

0.0400 (0.004849)

7

0.0410 (0.006454)

3

0.0140 (0.002410)

0.0950 (0.013713)

TD2

5

0.0381 (0.004699)

6

0.0353 (0.005079)

2

0.0096 (0.001313)

0.0830 (0.011091)

TD3

5

0.0420 (0.005231)

5

0.0331 (0.004745)

1

0.0049 (0.000631)

0.0800 (0.010607)

TD4

5

0.0420 (0.005231)

5

0.0331 (0.004745)

1

0.0049 (0.000631)

0.0800 (0.010607)

TD5

8

0.0590 (0.007020)

9

0.0578 (0.008782)

2

0.0094 (0.001426)

0.1262 (0.017228)

TD6

6

0.0490 (0.006118)

6

0.0350 (0.005175)

1

0.0022 (0.000540)

0.0862 (0.011833)

MTD1

7

0.0593 (0.006923)

9

0.0555 (0.008303)

3

0.0141 (0.002483)

0.1289 (0.017709)

MTD2

9

0.0709 (0.009138)

10

0.0620 (0.009262)

2

0.0104 (0.001590)

0.1433 (0.019990)

MTD3

9

0.0701 (0.008964)

10

0.0661 (0.009675)

2

0.0105 (0.001423)

0.1467 (0.020062)

MTD4

9

0.0701 (0.008965)

10

0.0661 (0.009674)

2

0.0104 (0.001423)

0.1466 (0.020062)

MTD5

10

0.0780 (0.010406)

11

0.0693 (0.010729)

2

0.0104 (0.001578)

0.1577 (0.022713)

MTD6

7

0.0561 (0.007525)

7

0.0464 (0.006590)

1

0.0034 (0.000642)

0.1059 (0.014757)

ACS Paragon Plus Environment

Page 17 of 31

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


Local energy densities (defined by Cremer and Kraka 65 ) at all critical points in the intermolecular region are reported in Table 4. They have positive values, as expected for non-covalent interactions. Larger value of electron density in a critical point of a certain type is characterized with larger energy density. Since localization of electron density between two monomers reflects interaction between the two units, it is to be expected that interaction energy could be represented as a function of ρ at intermolecular CPs. Certain empirical relations already exist. Zhikol et al. 26 analyzed stacked benzene dimers and presented stacking interaction energy as a function of the cosine of the angle formed by a cage critical point and ring critical points of aromatic rings, as well as electron density and its Laplacian in the CCP. Waller et al. 38 investigated a series of stacked aromatic species and presented interaction energy as dependent on the sum of electron densities in the critical points in the intermolecular region, since their results showed that higher values of electron density in the space between the two monomers corresponds to larger interaction energy. However, results presented in Table 4 show that for TD and MTD dimers increase in binding energy is not necessarily reflected in increase in total electron density in intermolecular CPs, nor in enlargement in local energy density in the corresponding points. Namely, previously analyzed dimers contained less atoms, and lower number of critical points: Zhikol et al. studied different conformations of benzene dimers, whereas Waller et al. analyzed monosubstituted benzene – benzene dimers with certain imposed constraints (for example, amino groups were kept planar, as they are in the case of a monomer). Obviously, full relaxation of dimeric forms results in appearance of additional critical points and uneven electron distribution, which does not lead to a simple relationship between binding energy and electron density (or local energy density) in the intermolecular region.

Ionized dimers In this section we will study influence of ionization on dimers’ stability and electron redistribution that accompanies it. Optimized structures of the ionized thymine dimers are given in the Supporting Information, Tables S13-S18. We used optimized geometries of neutral forms as initial 17 ACS Paragon Plus Environment


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

structures for optimizations. Ions’ structures do not significantly differ from the neutral ones, apart from TD5 which transformed into a T-shaped dimer with an N-H· · · O hydrogen bond (Supporting Information, Figure S1). In the following we will focus on changes that occur upon ionization of the most stable thymine dimer, TD1, and briefly discuss changes in other dimers. First, let us analyze influence of ionization on binding energies. Since BSSE correction is important for precise estimation of interaction energies in π -stacked systems, we first estimated charges of the monomers upon vertical ionization (TD1+ v is an ionized form of the structure that corresponds to the optimized neutral species) and within an optimized cationic species (TD1+ a ). In both cases, charges of both monomers amount to +0.5 e, which can be explained by the shape of the HOMO orbital, Figure 4, which is equally spread among both monomers both in TD1 and TD1+ a . As BSSE correction can be performed only provided charges of the fragments are integer numbers, we were not able to estimate it for the ionized species. Thus, we will only compare change in binding energy upon ionization (energy release when neutral thymine molecule T and cation T+ ) form a dimer by using energies of the dimers without counterpoise correction. This value amounts to -79 kJ/mol, which can be compared to -46 kJ/mol for the TD1 dimer without BSSE corrections. Significantly stronger interaction between the monomers in TD1+ a compared to TD1 reflects difference in charge distribution in neutral species and cations. A deeper insight into attractive interactions between monomers in TD1 and TD1+ a can be acquired from analysis of another property of electron density which proved to be significant for

TD1

TD1

+

Figure 4: HOMO orbital of the TD1 dimer (left) and the optimized cation TD1+ a (right).


Page 18 of 31

Page 19 of 31

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


Figure 5: Reduced density gradient isosurfaces for s=0.5 (ρ < 0.1 a.u.) with mapped function sign(λ2 ) · ρ in the range [-0.02 (blue), 0.02 (red)] a.u. Left to right: neutral T, neutral TD1, TD1+ v and TD1+ a. investigation of non-covalent interactions – reduced density gradient 39 (RDG) defined as

s=

1 2 3π 2

1/3

|∇ρ | . ρ 4/3

(2)

This dimensionless function is characterized with low values in regions with weak intramolecular interactions, although low values for ρ and s are specific not just for van der Waals interactions, but also for hydrogen bonding and steric repulsion. Identification of the type of weak interactions can be performed according to the sign of the middle eigenvalue of the Hessian of electron density, λ2 , which is negative for bonding interactions. Van der Waals interactions can be distinguished from hydrogen bonding by lower values of ρ in the bonding region. To incorporate both value of electron density and sign of λ2 , Yang and coworkers 39 proposed to map their product sign(λ2 ) · ρ on isosurfaces or RDG. Such a plot for T and TD1 is presented in Figure 5. Note that the chosen isosurface s=0.5 captures all intermolecular critical points (cf. Figure 3). Since deep blue/red color is used for more pronounced attractive/repulsive interactions that correspond to the mapping function equal to -0.02/0.02, green and light blue color demonstrate presence of weak attractive interactions (red color in the center of pyrimidine rings is due to steric repulsion). Also, small clouds are due to intramolecular interaction between neighboring carbonyl oxygen atoms and methyl groups, which are present in the monomer as well. Let us now discuss changes that occur upon vertical ionization (TD1+ v ), and subsequent relaxation (TD1+ a ). Although instantly after ionization takes place the charges of the monomer units increase by 0.5 e, as mentioned above, electron distribution in the intermolecular region does not undergo notable changes. However, relaxation on the ionic potential energy surface is accompa-



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

Page 20 of 31

nied with charge redistribution between two thymine species, which can clearly be seen in Figure 5 in the appearance of deep blue color. This is in accord with the larger value of binding energy for TD+ compared to TD discussed in the beginning of this section. The other five dimers display analogous changes of electron density distribution upon ionization. We will now turn to analysis of vertical and adiabatic ionization energies (IE). These values for thymine monomer amount to 9.25 and 8.95 eV, respectively. IEs for all six thymine dimers are compiled in Table 5. Vertical IEs in TDs are computed to be between 8.91 and 9.08 eV, i.e. lower than the value obtained for the monomer. Adiabatic IEs behave similarly. These results are in accord with values reported by Krylov and coworkers, 66,67 who analyzed ionization energies of selected nucleobase homodimers and heterodimers who concluded that both vertical and adiabatic energies are significantly affected by non-covalent interactions (stacking and hydrogen bonding). Concerning stacked dimers, the largest difference between IEs of the monomer and the dimers was noticed for the symmetric dimer. We arrived at the same conclusion, Table 5 – lowest adiabatic IE is observed for TD5, although vertical IE for this dimer is only 0.07 eV larger than for TD1 (and only 0.01 eV below that of TD6). + Table 5: Labels TD+ v and TDa stand for structures that correspond to vertical and adiabatic transitions. Vertical and adiabatic ionization energies (IE) are presented in eV, whereas sums of electron densities in critical points (∑ ρCP ) in the intermolecular region are shown in atomic units.

dimer

IE

TD1

∑ ρCP

dimer

0.095

TD4

IE

∑ ρCP 0.080

TD1+ v

8.91

0.087

TD4+ v

9.08

0.078

TD1+ a

8.60

0.114

TD4+ a

8.68

0.062

0.083

TD5

TD2

0.126

TD2+ v

9.05

0.083

TD5+ v

8.98

0.126

TD2+ a

8.68

0.140

TD5+ a

8.44

0.055

0.080

TD6

TD3

0.086

TD3+ v

9.08

0.078

TD6+ v

8.99

0.063

TD3+ a

8.68

0.062

TD6+ a

8.66

0.104


Page 21 of 31

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


Changes in the sum of electron densities in CPs in the intermolecular region are displayed in Table 5. Immediately after ionization, electron density does not significantly change in the intermolecular region. Depending on the nature of the dimer, relaxation of the ionized forms is marked either by increase (TD1, TD2, TD6) or decrease (TD3, TD4, TD5) of ∑ ρCP . This non uniform behavior also reflects complex nature of non-covalent interactions in thymine dimers.

Conclusions In this contribution we analyzed stability of stacked thymine dimers using density functional theory. It was shown that M06-2X, ω B97XD, B3LYP-D3 and CAM-B3LYP-D3 provide with comparable results, whereas PBE0-D3 did not manage to optimize one out of six structures. Comparison of binding energies in thymine dimers and 1-methylthymine dimers shows that methylation increases interaction energy. It was demonstrated that dimerization can take place at low temperatures, but not at room temperature due to considerable entropy decrease (T ∆S term). In each of six TDs we identified a low frequency mode that alters the distance between the two monomers and whose population might lead to dissociation. In one TD, this mode was not significantly coupled to other vibrational modes, whereas in other five TDs this low frequency mode was coupled to other vibrational degrees of freedom indicating complex vibrational dynamics. Born-Oppenheimer molecular dynamics simulations indicate that one of the analyzed dimers would not dissociate if prepared at 30 K and heated to room temperature. Bond, ring and cage critical points are identified in all investigated dimers. It was shown that increase in interaction energy is not reflected neither in larger number of critical points in the intermolecular region, nor in increase in total electron density (local energy density) computed at those critical points. Therefore, we did not succeed in designing a simple functional form that would enable computation of binding energy on the basis of total electron density in intermolecular critical points. This observation confirms complex interactions that lead to non-uniform distribution of electron density in the space between two monomer units.



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

Page 22 of 31

Our results show that ionization strengthens binding. This statement was confirmed by comparing binding energies of two neutral thymine molecules, and between a neutral and cation species, but also through analysis of reduced density gradient in the intermolecular region.

ASSOCIATED CONTENT Supporting Information M06-2X/6-311++G(2df,2pd) optimized geometries of stacked thymine dimers, 1-methylthymine dimers and cationic forms of thymine dimers, Figure of optimized TD5+ species. The Supporting Information is available free of charge on the ACS Publications website.

ACKNOWLEDGMENTS This work was financially supported by the Ministry of Education, Science, and Technological Development of the Republic of Serbia within the framework of the project 172040. Numerical simulations were run on the PARADOX cluster at the Scientific Computing Laboratory of the Institute of Physics Belgrade, supported in part by the Ministry of Education, Science, and Technological Development of the Republic of Serbia under project No. ON171017.

References (1) Kabeláˇc, M.; Hobza, P. Potential Energy and Free Energy Surfaces of All Ten Canonical and Methylated Nucleic Acid Base Pairs: Molecular Dynamics and Quantum Chemical ab Initio Studies. J. Phys. Chem. B 2001, 105, 5804–5817. (2) Cauët, E.; Rooman, M.; Wintjens, R.; Liévin, J.; Biot, C. Histidine-Aromatic Interactions in Proteins and Protein-Ligand Complexes: Quantum Chemical Study of X-ray and Model Structures. J. Chem. Theory Comput. 2005, 1, 472–483.


Page 23 of 31

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


(3) Rutledge, L. R.; Campbell-Verduyn, L. S.; Hunter, K. C.; Wetmore, S. D. Characterization of Nucleobase-Amino Acid Stacking Interactions Utilized by a DNA Repair Enzyme. J. Phys. Chem. B 2006, 110, 19652–19663. ˇ (4) Cerný, J.; Hobza, P. Non-covalent Interactions in Biomacromolecules. Phys. Chem. Chem. Phys. 2007, 9, 5291–5303. (5) Cysewski, P. A Post-SCF Complete Basis Set Study on the Recognition Patterns of Uracil and Cytosine by Aromatic and π -Aromatic Stacking Interactions with Amino Acid Residues. Phys. Chem. Chem. Phys. 2008, 10, 2636–2645. ˙ Zale´sny, R.; Czeleń, P. The Post-SCF Quantum Chemistry (6) Cysewski, P.; Czyz˙ nikowska, Z.; Characteristics of the Guanine-Guanine Stacking in B-DNA. Phys. Chem. Chem. Phys. 2008, 10, 2665–2672. (7) Churchill, C. D.; Wetmore, S. D. Noncovalent Interactions Involving Histidine: The Effect of Charge on π -π Stacking and T-Shaped Interactions with the DNA Nucleobases. J. Phys. Chem. B 2009, 113, 16046–16058. (8) Churchill, C. D.; Rutledge, L. R.; Wetmore, S. D. Effects of the Biological Backbone on Stacking Interactions at DNA-Protein Interfaces: The Interplay between the Backbone· · · π and π · · · π Components. Phys. Chem. Chem. Phys. 2010, 12, 14515–14526. (9) Parker, T. M.; Hohenstein, E. G.; Parrish, R. M.; Hud, N. V.; Sherrill, C. D. QuantumMechanical Analysis of the Energetic Contributions to π Stacking in Nucleic Acids versus Rise, Twist, and Slide. J. Am. Chem. Soc. 2013, 135, 1306–1316. (10) Hunter, R. S.; van Mourik, T. DNA Base Stacking: The Stacked Uracil/Uracil and Thymine/Thymine Minima. J. Comput. Chem. 2012, 33, 2161–2172. (11) Sherrill, C. D. Energy Component Analysis of π Interactions. Acc. Chem. Res. 2013, 46, 1020–1028. 23 ACS Paragon Plus Environment


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

(12) Riley, K. E.; Hobza, P. On the Importance and Origin of Aromatic Interactions in Chemistry and Biodisciplines. Acc. Chem. Res. 2013, 46, 927–936. (13) Song, Q.; Qiu, Z.; Wang, H.; Xia, Y.; Shen, J.; Zhang, Y. The Effect of Methylation on the Hydrogen-Bonding and Stacking Interaction of Nucleic Acid Bases. Struct. Chem. 2013, 24, 55–65. (14) Paytakov, G.; Gorb, L.; Stepanyugin, A.; Samiylenko, S.; Hovorun, D.; Leszczynski, J. Homodimers of Cytosine and 1-Methylcytosine. A DFT Study of Geometry, Relative Stability and H-NMR Shifts in Gas-Phase and Selected Solvents. J. Mol. Model. 2014, 20, 2115. (15) Angelina, E. L.; Andujar, S. A.; Tosso, R. D.; Enriz, R. D.; Peruchena, N. M. Non-Covalent Interactions in Receptor-Ligand Complexes. A Study Based on the Electron Charge Density. J. Phys. Org. Chem. 2014, 27, 128–134. (16) Estévez, L.; Sánchez-Lozano, M.; Mosquera, R. A. Understanding the Electron Density Reorganization upon Stacking vs. H-Bonding Interaction in Methyl Gallate-Caffeine Complexes. RSC. Adv. 2014, 4, 25018–25027. (17) Navarro, A.; Fernández-Liencres, M. P.; García, G.; Granadino-Roldán, J. M.; FernándezGómez, M. A DFT Approach to the Charge Transport Related Properties in Columnar Stacked π -Conjugated N-heterocycle Cores Including Electron Donor and Acceptor Units. Phys. Chem. Chem. Phys. 2015, 17, 605–618. (18) Parac, M.; Etinski, M.; Peric, M.; Grimme, S. A Theoretical Investigation of the Geometries and Binding Energies of Molecular Tweezer and Clip Host-Guest Systems. J. Chem. Theory Comput. 2005, 1, 1110–1118. (19) Sygula, A.; Fronczek, F. R.; Sygula, R.; Rabideau, P. W.; Olmstead, M. M. A Double Concave Hydrocarbon Buckycatcher. J. Am. Chem. Soc. 2007, 129, 3842–3843.


Page 24 of 31

Page 25 of 31

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


(20) Cho, Y.; Cho, W. J.; Youn, I. S.; Lee, G.; Singh, N. J.; Kim, K. S. Density Functional Theory Based Study of Molecular Interactions, Recognition, Engineering, and Quantum Transport in π Molecular Systems. Acc. Chem. Res. 2014, 47, 3321–3330. (21) Josa, D.; González-Veloso, I.; Rodríguez-Otero, J.; Cabaleiro-Lago, E. M. Tailoring Buckybowls for Fullerene Recognition. A Dispersion-Corrected DFT Study. Phys. Chem. Chem. Phys. 2015, 17, 6233–6241. (22) Gobre, V. V.; Tkatchenko, A. Scaling Laws for wan der Vaals Iterations in Nanostructured Materials. Nat. Commun. 2013, 4, 2341. (23) Müller-Dethlefs, K.; Hobza, P. Noncovalent Interactions: A Challenge for Experiment and Theory. Chem. Rev. 2000, 100, 143–167. (24) Lee, E. C.; Kim, D.; Jureˇcka, P.; Tarakeshwar, P.; Hobza, P.; Kim, K. S. Understanding of Assembly Phenomena by Aromatic-Aromatic Interactions: Benzene Dimer and the Substetuted Systems. J. Phys. Chem. A 2007, 111, 3446–3457. (25) Kim, K. S.; Karthikeyan, S.; Singh, N. J. How Different Are Aromatic π Interactions from Aliphatic π Interactions and Non-π Stacking Interactions. J. Chem. Theory Comput. 2011, 7, 3471–3477. (26) Zhikol, O. A.; Shishkin, O. V.; Lyssenko, K. A.; Leszczynski, J. Electron Density Distribution in Stacked Benzene Dimers: A New Approach Towards the Estimation of Stacking Interaction Energies. J. Chem. Phys. 2005, 122, 144104. (27) Li, H.; Lu, Y.; Liu, Y.; Zhu, X.; Liu, H.; Zhu, W. Interplay Between Halogen Bonds and π − π Stacking Interactions: CSD Search and Theoretical Study. Phys. Chem. Chem. Phys. 2012, 14, 9948–9955. (28) Gu, J.; Wang, J.; Leszczynski, J. Stacking and H-Bonding Patterns of dGpdC and dGpd-



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

CpdG: Performance of the M05-2X and M06-2X Minnesota Density Functionals for the Single Strand DNA. Chem. Phys. Lett. 2011, 512, 108–112. (29) Petković, M.; Ristić, M.; Etinski, M. Stability and Anharmonic N-H Stretching Frequencies of 1-Methylthymine Dimers: Hydrogen Bonding versus π -Stacking. J. Phys. Chem. A 2016, 120, 1536–1544. (30) Roca-Sanjuán, D.; Merchán, M.; Serrano-Andrés, L.; Rubio, M. Ab initio Determination of the Electron Affinities of DNA and RNA nucleobases. J. Chem. Phys. 2008, 129, 095104. (31) Bravaya, K. B.; Epifanovsky, E.; Krylov, A. I. Four Bases Score a Run: Ab Initio Calculations Quantify a Cooperative Effect of H-Bonding and π -Stacking on the Ionization Energy of Adenine in the AATT Tetramer. J. Phys. Chem. Lett. 2012, 3, 2726–2732. (32) Sun, H.; Zhang, S.; Sun, Z. Applicability of Optimal Functional Tuning in Density Functional Calculations of Ionization Potentials and Electron Affinities of Adenine-Thymine Nucleobase Pairs and Clusters. Phys. Chem. Chem. Phys. 2015, 17, 4337–4345. (33) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A Long-Range Correction Scheme for Generalized-Gradient-Approximation Exchange Functionals. J. Chem. Phys. 2001, 115, 3540–3544. (34) Sun, H.; Zhang, S.; Zhong, C.; Sun, Z. Theoretical Study of Excited States of DNA Base Dimers and Tetramers Using Optimally Tuned Range-Separated Density Functional Theory. J. Comput. Chem. 2016, 37, 684–693. (35) Hunter, C. A.; Sanders, J. K. M. The Nature of π -π Interactions. J. Am. Chem. Soc. 1990, 112, 5525–5534. (36) Cozzi, F.; Cinquini, M.; Annunziata, R.; Dwyer, T.; Siegel, J. S. Polar/π Interactions between Stacked Aryls in 1,8-Diarylnaphthalenes. J. Am. Chem. Soc. 1992, 114, 5729–5733.


Page 26 of 31

Page 27 of 31

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


(37) Paliwal, S.; Geib, S.; Wilcox, C. S. Molecular Torsion Balance for Weak Molecular Recognition Forces. Effects of "Tilted-T" Edge-to-Face Aromatic Interactions on Conformational Selection and Solid-State Structure. J. Am. Chem. Soc. 1994, 116, 4497–4498. (38) Waller, M. P.; Robertazzi, A.; Platts, J. A.; Hibbs, D. E.; Williams, P. A. Hybrid Density Functional Theory for π -Stacking Interactions: Application to Benzenes, Pyridines, and DNA Bases. J. Comput. Chem. 2006, 27, 491–504. (39) Johnson, E. R.; Keinan, S.; Mori-Sánchez, P.; Contreras-García, J.; Cohen, A. J.; Yang, W. Revealing Noncovalent Interactions. J. Am. Chem. Soc. 2010, 132, 6498–6506. (40) 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. et al., Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (41) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Account 2008, 120, 215–241. (42) McLean, A. D.; Chandler, G. S. Contracted Gaussian Basis Sets for Molecular Calculations. I. Second Row Atoms, Z=11-18. J. Chem. Phys. 1980, 72, 5639–5648. (43) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. XX. A Basis Set for Correlated Wave Functions. J. Chem. Phys. 1980, 72, 650–654. (44) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. v. R. Efficient Diffuse FunctionAugmented Basis-Sets for Anion Calculations. 3. The 3-21+G Basis Set for 1st-Row Elements, Li-F. J. Comp. Chem. 1983, 4, 294–301. (45) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-Consistent Molecular Orbital Methods. 25. Supplementary Functions for Gaussian Basis Sets. J. Chem. Phys. 1984, 80, 3265–3269. 27 ACS Paragon Plus Environment


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

(46) Gu, J.; Wang, J.; Leszczynski, J.; Xie, Y.; Schaefer III, H. To Stack or Not to Stack: Performance of a New Density Functional for the Uracil and Thymine Dimers. Chem. Phys. Lett. 2008, 459, 164–166. (47) van Mourik, T. Comment on ’To Stack or Not to Stack: Performance of a New Density Functional for the Uracil and Thymine Dimers’ [Chem. Phys. Lett. 459 (2008) 164]. Chem. Phys. Lett. 2009, 473, 206–208. (48) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615– 6620. (49) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652. (50) 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 1988, 37, 785–789. (51) Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51–57. (52) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (53) Perdew, J. P.; Burke, K.; Ernzerhof, M. Eratta: Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1997, 78, 1396. (54) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158–6169. (55) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures With Reduced Errors. Mol. Phys. 1970, 19, 553–566. 28 ACS Paragon Plus Environment

Page 28 of 31

Page 29 of 31

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


(56) Encyclopedia of Computational Chemistry A-D;

Schleyer, P., Allinger, N., Clark, T.,

Gasteiger, J., Kollman, P., Schaefer III, H., Schreiner, P., Eds.; John Wiley & Sons Ltd., West Sussex, 1998; p 13. (57) Petković, M. Infrared Spectroscopy of ClONO2 and BrONO2 Investigated by Means of Anharmonic Force Fields. Chem. Phys. 2007, 331, 438–446. (58) Keith, T. A. AIMAll, version 15.09.27, Professional; TK Gristmill Software, Overland Park KS, USA, 2014 (aim.tkgristmill.com). (59) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab initio Parametrization of Density Functional Dispersion Correction for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (60) Kim, K. S.; Mhin, B. J.; Choi, U.-S.; Lee, K. Ab initio Studies of the Water Dimer Using Large Basis Sets: The Structure and Thermodynamic Energies. J. Chem. Phys. 1992, 97, 6649–6662. (61) Tramer, A.; Jungen, C.; Lahmani, F. Energy Dissipation in Molecular Systems; Springer Berlin Heidelberg, 2005. (62) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press, Oxford, 1990. (63) The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design; Matta, C. F., Boyd, R. J., Eds.; WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim, 2007. (64) Bader, R. F. W. A Bond Path: A Universal Indicator of Bonded Interactions. J. Phys. Chem. A 1998, 102, 7314–7323. (65) Cremer, D.; Kraka, E. A Description of the Chemical Bond in Terms of Local Properties of Electron Density and Energy. Croat. Chem. Acta. 1984, 57, 1259–1281.



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

(66) Golubeva, A. A.; Krylov, A. I. The Effect of π -Stacking and H-Bonding on Ionization Energies of a Nucleobase: Uracil Dimer Cation. Phys. Chem. Chem. Phys. 2009, 11, 1303–1311. (67) Bravaya, K. B.; Kostko, O.; Ahmed, M.; Krylov, A. I. The Effect of π -Stacking, H-Bonding, and Electrostatic Interactions on the Ionization Energies of Nucleic Acid Bases: AdenineAdenine, Thymine-Thymine and Adenine-Thymine Dimers. Phys. Chem. Chem. Phys. 2010, 12, 2292–2307.


Page 30 of 31

Page 31 of 31

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


Graphical TOC Entry hn


Using Density Functional Theory To Study Neutral and Ionized

Recommend Documents