Table Of Content

First-principles GW calculations for DNA and RNA nucleobases Carina Faber1,2, Claudio Attaccalite1, Valerio Olevano1, Erich Runge2, Xavier Blase1 1Institut Néel, CNRS and Université Joseph Fourier, B.P. 166, 38042 Grenoble Cedex 09, France. 2Institut fu¨r Physik, Technische Universita¨t Ilmenau, 98693 Ilmenau, Germany. (Dated: January 20, 2011) Onthebasisoffirst-principlesGWcalculations, westudythequasiparticlepropertiesofthegua- nine, adenine, cytosine, thymine, and uracil DNA and RNA nucleobases. Beyond standard G0W0 calculations, starting from Kohn-Sham eigenstates obtained with (semi)local functionals, a simple 1 self-consistencyontheeigenvaluesallowstoobtainverticalionizationenergiesandelectronaffinities 1 within an average 0.11 eV and 0.18 eV error respectively as compared to state-of-the-art coupled- 0 clusterandmulti-configurational perturbativequantumchemistryapproaches. Further,GW calcu- 2 lationspredictthecorrectπ-characterofthehighestoccupiedstate,thankstoseverallevelcrossings betweendensityfunctionalandGWcalculations. Ourstudyisbasedonarecentgaussian-basis im- n plementation of GW with explicit treatment of dynamical screening through contour deformation a J techniques. 9 PACSnumbers: 31.15.A-,33.15.Ry,31.15.V- 1 ] h Thedeterminationoftheionizationenergies(IE),elec- brings the ionization energies in much better agreement p tronic affinities (EA) and character of the frontier or- withexperimentandhigh-levelquantumchemistrycalcu- - bitals of DNA and RNA nucleobases is an important lations. Theseresultsdemonstratetheimportanceofself- m step towards a better understanding of the electronic consistencyontheeigenvalueswhenperformingGWcal- e properties and reactivity of nucleotides and nucleosides culations in molecular systems starting from (semi)local h along the DNA/RNA chains. Important phenomena DFTfunctionals,andthemeritsofasimpleschemebased c . such as nucleobases/protein interactions, defining the on a G0W0 calculation starting from Hartree-Fock like s DNA functions1, or damages of the genetic material eigenvalues. c i through oxidation or ionizing radiations2, are strongly TheGWapproachisaGreen’sfunctionformalismusu- s related to these fundamental spectroscopic quantities. allyderivedwithinafunctionalderivativetreatment14,26 y h Even though nucleobases in DNA/RNA strands are allowing to prove that the two-body Green’s function p connected within the nucleotides to phosphate groups (G ), involved in the equation of motion of the one- 2 [ through a five-carbon sugar, several studies show that body time-ordered Green’s function G, can be recast the highest-occupied orbital (the HOMO level) in nu- intoanon-localandenergy-dependentself-energyopera- 1 v cleotides, which is responsible e.g. for the sensitivity of tor Σ(r,r′|E). This self-energy Σ accounts for exchange 8 themoleculetooxidationprocesses,remainslocalizedon and correlation in the present formalism. Since it is 3 the nucleobases3. Figure 1 shows the structures of the energy-dependent, it must be evaluated at the E = εQP i 7 DNA and RNA nucleobases, i.e. the purines - adenine quasiparticle energies, where (i) indexes the molecular 3 (A) andguanine (G), and the pyrimidines - cytosine (C) energy levels. This self-energy involves G(r,r′|ω), the 1. as well as thymine (T) in DNA and uracil (U) in RNA. dynamically-screenedCoulombpotentialW(r,r′|ω),and 0 the so-called vertex correction Γ. A set of exact self- Besides the overarching fundamental interest in un- 1 consistent (closed) equations connects G, W, Γ, and derstanding complex biological processes at the micro- :1 scopiclevel,ab initio calculationsofisolatednucleobases the independent-electron/full polarisabilities χ0(r,r′|ω) v and χ(r,r′|ω), respectively. In the GW approximation are interesting since recent high-level quantum chem- Xi istry calculations4–6 allow to rationalize the rather large (GWA), thethree-body vertexoperatorΓis setto unity, yielding the following expression for the self-energy: spread of experimental results concerning the electronic r a propertiesofthenucleobasesinthegasphase7–13,inpar- i ticular asdue to the existence ofseveralisomersfor gua- Σ(r,r′|E) = dωeiω0+G(r,r′|E+ω)W(r,r′|ω) 2π Z nine and cytosine6. Thus, these molecules offer a valu- afobrlmeamliseman14–t1o8efxoprliosroelatthede omregraintsicomf otlheecusloes-c,aallloedngGtWhe W˜(r,r′|ω) = Z dr1dr2v(r,r1)χ0(r1,r2|ω)W(r2,r′|ω), line of recent systematic studies of small molecules19 or χ (r,r′|ω) = (f −f )φ∗i(r)φj(r)φ∗j(r′)φi(r′) molecules such as fullerenes or porphyrins of interest for 0 i j ε −ε +ω±iδ X i j electronic or photovoltaic applications20–25. i,j In the present work, we study by means of first- wherev(r,r′)isthebare(unscreened)Coulombpotential principles GW calculations the quasiparticle properties andW˜ =W−v. The(ε ,φ )are“zeroth-order”one-body i i oftheDNAandRNAnucleobases,namelyguanine,ade- eigenstates. Followingthelargebulkofwork18devotedto nine,cytosine,thymineanduracil. Weshowinparticular GWcalculationsinsolids,surfaces,graphene,nanotubes, that the GW correction to the Kohn-Sham eigenvalues or nanowires, we use here Kohn-Sham DFT-LDA eigen- 2 FIG. 1: (Color online) Schematic representation of the molecular structure of (a) guanine (G9K), (b) adenine, (c) cytosine (C1), (d) thymine, and (e) uracil. Black, brown, red, white atoms are carbon, nitrogen, oxygen, and hydrogen, respectively. The G9K and C1 notations for theguanine and cytosine tautomers are consistent with Ref. 6. states. It is shown below, and in Refs. 19,25,27,28, that gaussianswith α=(0.1,0.4,1.5)a.u. describe hydrogen33. Hartree-Fock(orhybrid)solutionsmayconstitutebetter Ionization energies. We now comment on the val- starting points for molecular systems. (fi,fj) are Fermi- ues of the calculated first ionization energy (IE) as com- Dirac occupation numbers, and δ an infinitesimal such piled in the Table and Fig. 2. The comparison to the that the poles of W fall in the second and fourth quad- experimentaldataiscomplicatedbythe0.2-0.3eVrange rants of the complex plane. In the GW approximation, spannedbythevariousexperimentalreports(verticalar- the self-energy operator can be loosely interpreted as a rows Fig. 2). An additional complication in the case of generalization of the Hartree-Fock method by replacing cytosine and guanine, beyond the intrinsic difficulties in the bare Coulombpotentialwitha dynamicallyscreened accuratelymeasuringionizationenergiesinthegasphase, Coulomb interaction accounting both for exchange and is that several gas phase tautomers exist6 which differ (dynamical) correlations. An important feature of the fromthe so-calledC1-cytosineandG9K-guanineisomers GW approach is that not only ionization energies and commonlyfoundinDNA(seeFig.1). State-of-the-artab electronic affinities can be calculated, but also the full initio quantum chemistry calculations, namely coupled- quasiparticle spectrum. Further, both localized and in- cluster CCSD(T) and multiconfigurational perturbation finite systems can be treated on the same footing with (CASPT2)methods4,5,studiedthenucleobasetautomers long and short range screening automatically accounted that can be found along the DNA/RNA strands. More forintheconstructionofthescreenedCoulombpotential recently, equation of motion coupled-cluster techniques W. More details about the present implementation can (EOM-IP-CCSD) were performed on several isomers6. be found in Ref. 25. Allmethodsagreetowithin0.04eVfortheaverageIEof the A,G,C,T tautomersweconsiderhere,withamaxi- Ourcalculationsarebasedonarecentlydevelopedim- plementation of the GW formalism (the Fiesta code) mumdiscrepancyof0.09eVin the case ofthymine. The CASPT2andCCSD(T)calculationsagreetowithin0.03 using a gaussianauxiliary basis to expand the two-point operators such as the Coulomb potential, the suscepti- eV for all molecules. These theoreticalIE are commonly bilities or the self-energy25. Dynamical correlations are consideredasthemostreliablereferencesandlandwithin included explicitly through contour deformation tech- the experimentalerrorbars,exceptfor the cytosine(C1) niques. We start with a ground-state DFT calculation case where the calculated IEs are slightly smaller than using the Siesta package29 and a large triple-zeta with the experimental lower bound34 (see Table and Fig. 2). double polarization (TZDP) basis30. We fit the radial Clearly, the ionization energy within DFT-LDA, as partofthenumericalbasisgeneratedbytheSiestacode given by the negative HOMO Kohn-Sham level energy, by up to five contracted gaussians in order to facilitate significantly underestimates the IE by an average of the calculation of the Coulomb matrix elements and of ∼2.5 eV (29%)35. The self-energy correction at the the matrix elements hφ |β|φ i of the auxiliary basis (β) G W (LDA) level improves very significantly the situa- i j 0 0 between Kohn-Shamstates. Such a scheme allowsto ex- tionbybringingthe errorto anaverage0.5eV(5.7%)as ploit the analytic relations for the products of gaussian compared to state-of-the-art quantum chemistry results. orbitals centered on different atoms or for their Fourier However, as emphasized in recent papers19,25,27,28, the transform25. Our auxiliary basis for first row elements overscreeninginducedbystartingwithLDAeigenvalues, is the tempered basis31 developed by Kaczmarski and which dramatically underestimate the band gap, tends coworkers32. Suchabasiswastestedrecentlyinasystem- to produce too small ionization energies. This problem aticstudyofseveralmoleculesofinterestforphotovoltaic canbe solvedatleastpartlybyperforminga simpleself- applications25. Fourgaussiansforeachl-channelwithlo- consistency on the eigenvalues. We shall refer to this calization coefficients α=(0.2,0.5,1.25,3.2) a.u. are used approach as GW henceforth. Such a self-consistency on forthe(s,p,d)channelsofC,O,andNatoms,whilethree the eigenvalues leads to a much reduced average error 3 10 LDA QuantChem not quite clear what should be the best starting point. 9.5 G0 W0 (LDA) Exp. range Next, we address the character of the HOMO level GW of cytosine and uracil. It changes from DFT-LDA to GW calculations. We plot in Fig. 3(a-d) the C1- 9 ) cytosine DFT-LDA Kohn-Sham HOMO to (HOMO-3) V (e8.5 eigenstates. TheLDAHOMOlevelisaninplaneσ state gy witha strongcomponentonthe (px,py) oxygenorbitals. er 8 Suchastateis labeledσ inthe Tableandinthefollow- n O Ionization e67..755 enin enined enisotyC enimyhT licarU wizGσneOeWgni.geshTtrataihnntoedegn(iπHGtshO0mpeWMuoosl0xOeh(yceH-gud1Fel)andtlreio(avopgear)zlb)isaistiopgaarpnlb.mriifiotWoacaralciehntahstenltsiydna,natltdohhdaweeeredlGLroπDc0ea-WnAslietz0areH(gtLdeyODwbMaAeintnO)hd-, au A the π state becomes the HOMO level. This level cross- G 6 ingbringstheGWcalculationsinagreementwithmany- body quantum chemistry calculations, which all predict 5.5 the π state to be the HOMO level. The same levelcross- ingisobservedinthecaseofuracilwiththeLDAHOMO FIG. 2: (Color online) Ionization energies in eV. The ver- and(HOMO-1)levelsbeingσ andπ-statesrespectively, tical (maroon) error bars indicate the experimental range. O whileallGWresultsandquantumchemistrycalculations Triangles up (light blue): LDA values; (green) squares: G0W0(LDA) values; full black diamond: GW values; (red) predictareverseordering. Ourinterpretationisthatthe empty circles (QuantChem abbreviation): quantum chem- very localized σO state suffers much more from the spu- istry,namelyCCSD(T),CASPT2andEOM-IP-CCSD,values rious LDA self interaction than the rather delocalized π (see text). state. Even though it would be wrong to reduce the dy- namicalGWself-energyoperatortoaself-interactionfree functional,theGWcorrectioncertainlycuresinpartthis well-known problem. The other bases, namely guanine, of 0.11 eV (∼1.3%) as compared to the quantum chem- adenine, and thymine, all show the correct π-character istry reference. This good agreement certainly indicates for the HOMO level. the reliability of the present GW scheme for such sys- The HOMO to (HOMO-1) energy difference averages tems. As shown in Fig. 2, the largest discrepancies are to 0.80 eV and 1.12 eV within CASPT2 and EOM-IP- observedforguanineandadenine(thepurines),whilethe CCSD, respectively. Clearly, the average LDA energy agreement is excellent for the three remaining bases. In recent work, it was shown that for small molecules a non-self-consistent G W calculation starting from 0 0 (a) (b) Hartree-Fock eigenstates leads for the ionization energy to better results than a full self-consistent GW calculation where the wavefunctions are updated as well19,27. Consistentwiththisobservation,asimpleschemerelying onanHartree-Fock-likeapproachwassuccessfullytested onsilane,disilane,andwater28,andlargermoleculessuch asfullerenesorporphyrins25. Inthis“G W onHartree- 0 0 Fock (HF)” ansatz, the input eigenvalues (ǫ˜ ) are com- (σ ) (π) n O puted within a diagonal first-order perturbation theory where the DFT exchange-correlationcontribution to the (c) (d) eigenvalues is replaced by the Fock exchange integral, namely: ǫ˜ =ǫLDA+<ψLDA|Σ −VLDA|ψLDA >. n n n x xc n where Σ is the Fock operator. This approach, labeled x G W (HF ) in the Table, produces an average error (σ) (π´) 0 0 diag of 0.22 eV (∼2.6%). This good agreementwith both the GW and quantum chemistry calculations clearly speaks FIG. 3: (Color online) Isodensity surface plot of the HOMO in favor of this simple scheme for molecular systems, or (σO), HOMO-1 (π), HOMO-2 (σ), and HOMO-3 (π′) LDA the full G W (HF) calculations tested in Ref. 19, which Kohn-Shameigenstatesofcytosine. WithinGW,theordering 0 0 alsoavoidsseekingself-consistency. Adifficultissuelying of states becomes π,π′,σO, σ for HOMO to HOMO-3 (see aheadconcernse.g. hybridsystems,suchassemiconduct- text). ing surfaces grafted by organic molecules, for which it is 4 Vertical ionization energies and vertical electronic affinities LDA-KS G0W0(LDA) GW G0W0(HFdiag) CASa,b/CCa,b EOMc Experimentd,e,f,g G-LUMO 1.80 -1.04 -1.58 -1.77 -1.14a/ G-HOMO 5.69 7.49 7.81 7.76 8.09b/8.09b 8.15 8.0-8.3d G-HOMO-1 6.34 8.78 9.82 9.78 9.56b/ 9.86 9.90g A-LUMO 2.22 -0.64 -1.14 -1.30 -0.91a/ -0.56 to -0.45e A-HOMO 6.02 7.90 8.22 8.23 8.37b/8.40b 8.37 8.3-8.5d, 8.47f A-HOMO-1 6.28 8.75 9.47 9.51 9.05b/ 9.37 9.45f C-LUMO 2.57 -0.45 -0.91 -1.05 -0.69a/-0.79a -0.55 to -0.32e C-HOMO 6.167 (σO) 8.21 (π) 8.73 (π) 9.05 (π) 8.73b (π)/8.76b 8.78 (π) 8.8-9.0d, 8.89f C-HOMO-1 6.172 (π) 8.80 (σO) 9.52 (π’) 9.87 (π’) 9.42b (σO)/ 9.54 (π’) 9.45g, 9.55f C-HOMO-2 6.806 (σ) 8.92 (π’) 9.89 (σO) 10.36 (σO) 9.49b (π’)/ 9.65 (σO) 9.89f C-HOMO-3 6.809 (π′) 9.38 (σ) 10.22 (σ) 10.64 (σ) 9.88b(σ)/ 10.06 (σ) 11.20f T-LUMO 2.83 -0.14 -0.67 -0.77 -0.60a/-0.65a -0.53 to -0.29e T-HOMO 6.54 8.64 9.05 9.05 9.07b/9.04b 9.13 9.0-9.2d, 9.19f T-HOMO-1 6.68 9.34 10.41 10.40 9.81b/ 10.13 9.95-10.05d,10.14f U-LUMO 3.01 -0.11 -0.64 -0.71 -0.61a/-0.64a -0.30 to -0.22e U-HOMO 6.72 (σO) 9.03 (π) 9.47 (π) 9.73 (π) 9.42b (π)/9.43b 9.4-9.6d U-HOMO-1 6.88 (π) 9.45 (σO) 10.54 (σO) 10.96 (σO) 9.83b (σO)/ 10.02-10.13d U-HOMO-2 7.55 (σ) 9.88 (π’) 10.66 (π’) 11.06 (π’) 10.41b (π’)/ 10.51-10.56d U-HOMO-3 7.66 (π’) 10.33 (σ) 11.48 (σ) 11.90 (σ) 10.86b (σ)/ 10.90-11.16d MAELUMO 3.29 0.33 0.18 0.31 MAEHOMO 2.5 0.5 0.11 0.22 TABLE I: Vertical ionization energies and electronic affinities in eV as obtained from the negative Kohn-Sham eigenvalues (LDA-KS),from non-self-consistent G0W0(LDA)calculations, from a GW calculation with self-consistency on theeigenvalues (GW), and from a non-self-consistent G0W0(HFdiag) calculation starting from Hartree-Fock-like eigenvalues. The σ or π character of the wavefunctions is indicated when the GW correction changes the level ordering as compared to DFT-LDA (seetext). TheacronymsCAS,CCandEOMstandforCASPT2,CCSD(T) andequationofmotion coupled-clusterhigh-level manybodyquantumchemistrycalculations,respectively. TheoreticalvaluesarereportedfortheC1-cytosineandG9K-guanine, while theexperimental values average over several tautomers. The MAE is themean absolute error in eV as compared to the quantum chemistry reference calculations in columns 6 and 7. aRef. 5. bRef. 4. cRef. 6. dCompiled in Ref. 4. eCompiled in Ref. 5. fRef. 10. gRef. 8. spacing of 0.22 eV is significantly too small. We find GWelectronicaffinitiesarequitesatisfying,withaMAE that the 0.77 eV G W (LDA) average value is close to of 0.18 eV. Such an agreement is rather impressive since 0 0 the CASPT2 results, while the larger 1.29 eV GW re- the LDA electronic affinities show the wrong sign, with sult falls closer to the EOM-IP-CCSD energy difference. a discrepancyas comparedto CASPT2 rangingfrom2.9 Averaging over all isomers, the experimental HOMO eV to 3.6 eV. We observe that while the G W EAs are 0 0 to (HOMO-1) energy spacing comes to 0.97 eV, in be- smaller (in absolute value) than the quantum chemistry tween the G W (LDA) or CASPT2 results and the GW ones, the GW EAs are larger. This contrasts with the 0 0 or EOM-IP-CCSD values. Even though it is too early IE case where both G W and GW values were smaller 0 0 for final conclusions about the merits of the various ap- (see Fig. 2). Similar to the quantum chemistry case, the proaches,itseemsfairtostatethattheLDAvalueissig- GW values are found to systematically overestimate the nificantly too small, and that the situation is improved experimental results. Further study is needed to under- significantly by the GW correction. stand such a discrepancy between theoretical and avail- Electronic affinities. We conclude this study by ex- able experimental results. ploring the electronic affinity (EA) of the nucleobases. Inconclusion,wehavestudiedonthe basisofab initio TheyareprovidedintheTableasthenegativesignofthe GW calculations the ionization energies and electronic LUMOKohn-Shamenergies. Experimentaldataforgua- affinities of the DNA and RNA nucleobases, guanine, nine are missing. Further, the CASPT2 and CCSD(T) adenine, cytosine, thymine and uracil. While a stan- results5 are clearly larger (in absolute value) than the dard G W (LDA) calculation yields ionization energies 0 0 highest experimental estimates. While again part of the that are 0.5 eV awayfromCCSD(T)/CASPT2 reference discrepancy may come from the presence of several tau- quantum chemistry calculations, self-consistency on the tomers in the gas phase, it certainly results as well from eigenvalues brings the agreementto an excellent 0.11 eV thefactthattheelectronicaffinityisnegative. Adetailed averageabsoluteerror. AsimpleG W calculationstart- 0 0 discussiononthe experimentaldifficulties in probingun- ingfromHartree-Fock-likeeigenvalues,avoidingtheneed bound states is presented in Ref. 6. Taking again the for self-consistency,shifts the agreementto 0.22eV. The CCSD(T) and CASPT2 calculations5 as a reference, the possibilityofbringingthecalculatedvaluestowithin0.1- 5 0.2eVfromstate-of-the-artreferencecalculationswitha biological systems in general. scheme, the GW formalism, which allows to treat both Acknowledgements. C.F. is indebted to the Euro- finite size and extended systems with a N4 scaling, and pean Union Erasmus program for funding. Calculations permits to obtain the full quasiparticle spectrum, paves havebeenperformedontheCIMENTplatforminGreno- the way to further studies of larger DNA strands and ble thanks to the Nanostar RTRA project. 1 Awell-knownexampleisthatofthe“replication”proteins al.,Phys. Chem. Chem. Phys. 12, 10817 (2010). thatcopythenucleobasesinDNAtranscriptionandrepli- 25 X.Blase, C.AttaccaliteandV.Olevano,arXiv:1011.3933. cation. See e.g. A. Travers, in DNA-Protein Interactions. 26 P. C. Martin and J. Schwinger, Phys. Rev. 115, 1342 Springer (1993); P.B. Dervan,Science 232. 464 (1983). (1959). 2 H. Lodish, A. Berk, P. Matsudaira, C. A. Kaiser et al., in 27 K. Kaasbjerg, K. S. Thygesen, Phys. Rev. B 81, 085102 Molecular Biology of the Cell, WH Freeman: New York, (2010). NY.5th ed. (2004). 28 P. H.Hahn,W. G. Schmidtand F. Bechstedt, Phys. Rev. 3 D. M. Close and K. T. Øhman, J. Phys. Chem. A 112, B 72, 24545 (2005). 11207 (2008); and references therein. 29 J.M.Soleret al.,J.Phys.: Condens.Mater14,2745-2779 4 D. Roca-Sanjuan, M. Rubio, M. Merchan et al., J. Chem. (2002). Phys. 125, 084302 (2006). 30 It was shown in Refs. 19,25 for a large set of molecules 5 D.Roca-Sanjuan,M.Merchan,L.Serrano-Andres,M.Ru- thatastandarddouble-zetapluspolarization (DZP)basis bio, J. Chem. Phys. 129, 095104 (2008). to expand the Kohn-Sham eigenstates yields already IE 6 K. B. Bravaya, O. Kostko, S. Dolgikh, A. Landau et al., within ∼ 0.1 eV as compared to the much larger TZDP J.Phys.Chem.A114,12305-12317(2010);andreferences basis.ThiscertainlyindicatesthatthepresentTZDPbasis therein. is an extremely good basis for such calculations. 7 N.S.Hush,A.S.Cheung,Chem.Phys.Lett.34,11(1975). 31 Forarecentanalysis, see: I.Cherkes,S.Klaiman,N.Mis- 8 D. Dougherty, E. S. Younathan, R. Voll et al., J. Elec. eyev, Int. J. Quant. Chem. 109, 2996 (2009); and refer- Spectro. Relat. Phenomena 13, 379 (1978). ences therein. 9 K.W.Choi,J.H.Lee,S.K.Kim,J.Am.Chem.Soc.127, 32 M. S. Kaczmarski, Y. C. Ma and M. Rohlfing,Phys. Rev. 15674 (2005). B 81, 115433 (2010). 10 A.B.Trofimov,J.Schirmer,V.B.Kobychev,A.W.Potts 33 We have tested the use of a larger and more diffuse even- et al., J. Phys. B: At.Mol. Opt.Phys. 39, 305 (2006). temperedbasiswithdecaycoefficientsα=(0.15,0.32,0.69, 11 M. Schwell, H. W. Jochims, H. Baumgartel et al., Chem. 1.48,3.2)a.u.forC,O,andNatoms.WefindthattheGW Phys. 353, 145 (2008). ionization energies and electronic affinities change by 0.08 12 I.L.Zaytseva,A.B.Trofimov,J.Schirmeret al.,J.Phys. eV and 0.015 eV (MAE) respectively with no systematic Chem. A 113, 15142 (2009). trend. 13 O. Kostko, K. Bravaya, A. Krylov, et al., Phys. Chem. 34 TheEOM-IP-CCSDcalculationspredictthattheIPofthe Chem. Phys. 12, 2860 (2010); and references therein. most stable guanine (G7K) and cytosine (C2b) isomers 14 L. Hedin,Phys. Rev.139, A796 (1965). (the most abundant in the gas phase) are 0.14 eV and 15 G. Strinati, H. J. Mattausch, W. Hanke, Phys. Rev. Lett. 0.08 eV larger than the IP of the G9K and C1 isomers 45, 290 (1980); idem,Phys. Rev.B 25, 2867 (1982). we study. Our calculations yields differences of 0.08 eV 16 M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 and 0.04 eV respectively (GW value). For cytosine, our (1986). GW calculations predict that the largest IE is that of the 17 R. W. Godby,M. Schlu¨ter, and L. J. Sham, Phys. Rev.B (C3a) tautomer which lands 0.19 eV above that of (C1). 37, 10159 (1988). Thisisconsistentwiththe0.2eVexperimentalrange,even 18 G. Onida,L.Reining, A.Rubio,Rev.Mod. Phys.74, 601 thoughlargerthanthe0.12eVdifferencereportedinRef.6 (2002). (EOM-IP-CCSD calculations). Clearly, averaging over all 19 C.Rostgaard,K.W.Jacobsen,K.S.Thygesen,Phys.Rev. isomers would bring our results in better agreement with B 81, 085103 (2010). experiment. We note however that such an average would 20 N. Dori et al.,Phys. Rev.B 73, 195208 (2006). require to know the abundance of each tautomer at the 21 M. L. Tiago, P. R. C. Kent, R.Q. Hood, F. A. Reboredo, experimentaltemperature.Thisstandsbeyondthepurpose J. Chem. Phys.129, 084311 (2008). of the present paper. 22 P. Umari, G. Stenuit,S. Baroni, Phys. Rev.B 79, 201104 35 VerysimilarresultsareobtainedwiththePBEfunctional. (R) (2009). J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 23 M. Palummo et al., J. Chem. Phys. 131, 084102 (2009). 77, 3865 (1996). 24 G. Stenuit, C. Castellarin-Cudia, O. Plekan, V. Feyer et