Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Electronic structure and spectroscopy in the gas and condensed phase: methodology and applications
(USC Thesis Other)
Electronic structure and spectroscopy in the gas and condensed phase: methodology and applications
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
ELECTRONICSTRUCTUREANDSPECTROSCOPYINTHEGASAND CONDENSEDPHASE:METHODOLOGYANDAPPLICATIONS. by VitaliiVanovschi ADissertationPresentedtothe FACULTYOFTHEGRADUATESCHOOL UNIVERSITYOFSOUTHERNCALIFORNIA InPartialFulfillmentofthe RequirementsfortheDegree DOCTOROFPHILOSOPHY (CHEMISTRY) May2009 Copyright 2009 VitaliiVanovschi Acknowledgements I’d like to thank my scientific adviser and mentor Prof. Anna I. Krylov. Her advice andsupportwerethecriticalsuccessfactorsintheprojectspresentedinthisthesis. Her enthusiasm toward the results I’ve obtained was highly encouraging. I appreciate the mentoring on the number various topics Anna provided me with. Also, I would like to thankAnnafortheteam-buildingactivitiessheorganized: hikes,trips,rock-climbing. Thanks to the current and the former members of Anna’s group, especially Vadim Mozhayskiy, Dr. Kadir Diri, Evgeny Epifanovsky, Dr. Sergey V. Levchenko, Lucasz KoziolandDr. PiotrA.Pieniazekfortheirhelpthroughttheprojects. I’dliketothankProf. CurtWittingforthefascinatingphysicaltheoriesI’velearned inhisclassandformentoringmeduringthecourseofstudies. I am grateful to our main collaborator on the effective fragment potential project, Prof. Lyudmila Slipchenko for answering all the numerous questions I had about this complicated method. Also I’d like to thank Prof. Mark S. Gordon for organizing a welcomevisittoIowaState Universitywhere, workinginthelaboratoryoftheoriginal EFPinventors,Imadeasubstantialprogressonthisproject. I’d like to thank Dr. Yihan Shao and Dr. Ryan Steele who helped me to figure out thenutsandboltsofQ-Chem. EspeciallyI’dliketothankDr. ChristopherWilliamsfor implementing the multipole integrals, required for the EFP electrostatics and polariza- tioncomponents. ii TableofContents Acknowledgements ii ListofFigures v ListofTables viii Abstract ix Chapter1: ImplementationoftheEffectiveFragmentPotentialmethod 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Electrostatics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3 Polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4 Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.5 Exchange-repulsion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.6 ImplementationofEFPinQ-Chem . . . . . . . . . . . . . . . . . . . . 28 1.7 EFPapplication: Beyondbenzenedimer . . . . . . . . . . . . . . . . . 37 1.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 Chapter1References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Chapter2: Structure,vibrationalfrequencies,ionizationenergies,andphoto- electronspectrumofthepara-benzyneradicalanion 57 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.2 Molecular orbital picture and the origin of vibronic interactions in the para-benzyneanion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 2.3 Computationaldetails . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 2.4 Equilibrium structures, charge distributions, and the D 2h → C 2v poten- tialenergyprofilesatdifferentlevelsoftheory . . . . . . . . . . . . . . 67 2.5 Vibrationalfrequenciesandthephotoelectronspectrumofpara-C 6 H − 4 . 77 2.6 DFTself-interactionerrorandvibronicinteractions . . . . . . . . . . . 82 2.7 Calculationofelectronaffinitiesofp-benzyne . . . . . . . . . . . . . . 86 2.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 Chapter2References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 iii Chapter3: Futurework 96 3.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 3.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Chapter3References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 Bibliography 105 Appendices 118 AppendixA: GAMESSinputforgeneratingEFPparametersofbenzene . . 118 AppendixB: Geometriesofbenezeneoligomers . . . . . . . . . . . . . . . 120 AppendixC: ConverteroftheEFPparameters . . . . . . . . . . . . . . . . 123 AppendixD: Fragment-fragmentelectrostaticinteraction . . . . . . . . . . 129 AppendixE: Buckinghammultipoles . . . . . . . . . . . . . . . . . . . . . 131 AppendixF: TheQ-CheminputforEFPwatercluster . . . . . . . . . . . . 135 AppendixG: Q-Chem’skeywordsaffectingtheEFPparameterscomputation 139 iv ListofFigures 1.1 Comparisonofcontinuum(a)anddiscrete(b)modelsforawatermolecule intheenvironmentofotherwatermolecules. . . . . . . . . . . . . . . 3 1.2 Uracil in aqueous solution. Active region (uracil) is shown with a balls andsticksmodel. Effectivefragments(water)aredepictedbylines. . . 4 1.3 Distributedmultipoleexpansionpointsforawatermolecule. . . . . . . 10 1.4 Sketchofthechargepenetrationphenomenonfortwowatermolecules. 13 1.5 The shape of the grid (monotonous gray area) for fitting the screening parametersofawatermolecule. . . . . . . . . . . . . . . . . . . . . . 15 1.6 Distributedpolarizabilitiesforawatermolecule. . . . . . . . . . . . . 16 1.7 Computing the induced dipoles through the two level self-consistent iterative procedure. Here Ψ is the ab initio wavefunction, μ is a set ofinduceddipoles,F Ψ istheforceduetotheabinitiowavefunctionand F μ istheforceduetotheinduceddipoles. . . . . . . . . . . . . . . . . 19 1.8 Distributeddispersionpointsforawatermolecule. . . . . . . . . . . . 23 1.9 Theexchange-repulsionLMOcentroidsforawatermolecule. . . . . . 27 1.10 ClassdiagramfortheEFPenergypartoftheimplementation. . . . . . 29 1.11 TheinternalstructureoftheFragmentParamsclass. . . . . . . . . . . 29 1.12 RepresentingtheorientationwiththeEulerangles. xyzisthefixedsys- tem,XYZistherotatedsystem,Nistheintersectionbetweenthexyand XYplanescalledthelineofnodes. Transformationbetweenthexyzand XYZsystemsisgivenbyαrotationinthexyplane,followedbyβ rota- tioninthenewy’z’plane(aroundtheNaxis),followedbyγ rotationin thenewx”y”plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 1.13 TherotationmatrixcorrespondingtotheEuleranglesα,β,γ. . . . . . 31 v 1.14 TheflowofcontrolbetweenQ-ChemandEFPenergymodule. . . . . . 32 1.15 Theformatforthe$efp paramssection. . . . . . . . . . . . . . . . . . 35 1.16 Theformatforthe$efp fragmentssection. . . . . . . . . . . . . . . . 35 1.17 Sandwitch (S), t-shaped (T) and parallel-displaced (PD) structures of thebenzenedimer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 1.18 Thestructuresofeightbenzenetrimers. . . . . . . . . . . . . . . . . . 45 1.19 Dependence of the non-additive energy share on the length of benzene oligomers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.1 Frontier MOs and leading electronic configurations of the two lowest electronicstatesatD 2h andC 2v structures. . . . . . . . . . . . . . . . . 62 2.2 ThetwolowestelectronicstatesofC 6 H − 4 canbedescribedby: (i)EOM- IP with the dianion reference; (ii) EOM-EA with the triplet neutral p- benzyne reference; (iii) EOM-SF with the quartet p-benzyne anion ref- erence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.3 D 2h structuresoptimizedatdifferentlevelsoftheory. . . . . . . . . . . 68 2.4 CCCangles,bondslengthscalculatedbydifferentcoupled-clustermeth- odsandB3LYP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 2.5 The distance between the radical-anion centers calculated by different coupled-clustermethodsandB3LYP. . . . . . . . . . . . . . . . . . . . 70 2.6 The NBO charge at C 1 (top panel) and dipole moments (bottom panel) at the optimized C 2v geometries. The ROHF-CCSD and ROHF-EA- CCSD values are calculated at the geometries optimized by the corre- spondingUHF-basedmethod. . . . . . . . . . . . . . . . . . . . . . . 71 2.7 Energy differences between the C 2v and D 2h structures (6-311+G** basis). Upperpanel: calculatedusingtheUHF-CCSDoptimizedgeome- tries. MCSCF and CASPT2N employ the MCSCF geometries. Lower panel: calculatedattherespectiveoptimizedgeometries. . . . . . . . . 73 2.8 ThepotentialenergyprofilesalongtheD 2h - C 2v scans. AttheD 2h struc- ture,CCSDandCCSD(T)energiesarecalculatedusingsymmetry-broken andpure-symmetryUHFreferences(seesection2.3). . . . . . . . . . . 75 2.9 Displacementsalongthenormalmodesuponionization(D 2h structure). 80 vi 2.10 Top: Singlet (solid line) and triplet (dotted line) bands of the electron photodetachmentspectrumforD 2h usingtheCCSDequilibriumgeome- triesandfrequenciesoftheanion. Bottom: theexperimental(solidline) andthecalculated(dottedline)spectraforD 2h usingtheCCSDstructure andnormalmodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 2.11 Theexperimental(solidline)andthecalculated(dottedline)spectrafor C 2v usingtheCCSDstructureandnormalmodes. . . . . . . . . . . . . 83 2.12 Theexperimental(solidline)andthecalculated(dottedline)spectrafor D 2h usingtheB3LYPstructureandnormalmodes. . . . . . . . . . . . 84 2.13 Energy differences between the D 2h and C 2v minima (using the CCSD structures)calculatedusingdifferentdensityfunctionals. . . . . . . . . 87 2.14 Adiabaticelectronaffinitiescalculatedbyenergydifferencesforthesin- glet (top panel) and triplet (bottom panel) state of the para-benzyne diradical. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 vii ListofTables 1.1 Thegeometricparametersandenergiesofthebenzenedimers . . . . . 41 1.2 Dissociationenergiesofthebenzenedimersatdifferentlevelsoftheory. 42 1.3 Total and many-body interaction energies of the benzene trimers (first group). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 1.4 Total and many-body interaction energies of the benzene trimers (sec- ondgroup). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 1.5 BindingenergiesofbenzeneoligomerscomputedwithEFP. . . . . . . 49 2.1 Energy differences between the C 2v and D 2h minima (eV) in the 6- 311+G** basis set. Negative values correspond to the C 2v minimum beinglower. ZPEsarenotincluded. . . . . . . . . . . . . . . . . . . . 74 2.2 Vibrationalfrequenciesofthepara-benzyneradicalanionandthedirad- ical. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 2.3 Displacementsinnormalmodesuponionization,A∗ √ amu. . . . . . . 80 2.4 Electronaffinities(eV)ofthea 3 B 1u andX 1 A g statesofthepara-benzyne diradical. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 viii Abstract Methoddevelopmentsandapplicationsofthecondensedphaseandgasphasemodeling techniquesarepresentedintwopartsofthiswork. Inthefirstpart,thetheoryoftheeffectivefragmentpotential(EFP)methoddesigned foraccuratemodelingofthemolecularpropertiesoftheextendedsystemsisintroduced. TheimplementationofEFPdevelopedinthisworkwithintheQ-Chemelectronicstruc- ture package is discussed in details. The code is applied to the studies of π−π inter- actions in benzene oligomers. The primary interest of these studies are the many-body and total non-additive π−π interactions. A detailed comparison of EFP with ab initio theories is presented for the set of eight benzene trimers. Three-body intermolecular interactionenergiesattheEFPleveloftheoryarewithin0.04kcal/molfromthefull ab initio results. To elucidate the asymptomatic behavior of the total non-additive binding energy, three types of linear benzene oligomers with 10, 20, 30 and 40 monomers are modeledusingournewEFPimplementation. Theresultofthisinvestigationisthatwith the increase in the size of linear benzene clusters, the share of the total non-binding energy approaches a constant and does not exceed 4% for the structures considered in thestudy. Thisresultindicatesthattherearenosubstantialnon-localeffectsinthebind- ingmechanismoflinearbenzeneclusters. In the second part, the equilibrium structure, vibrational frequencies, and ioniza- tion energies of the para-benzyne radical anion are characterized by coupled-cluster ix andequation-of-motionmethods. Vibronicinteractionswiththelow-lyingexcitedstate result in a flat potential energy surface along the coupling mode and even in lower- symmetryC 2v structures. AdditionalcomplicationsariseduetotheHartree-Fockinsta- bilitiesandnear-instabilities. Themagnitudeofvibronicinteractionswascharacterized bygeometricalparameters,chargelocalizationpatternsandenergydifferencesbetween theD 2h andC 2v structures. TheobservedtrendssuggestthattheC 2v minimumpredicted byseveraltheoreticalmethodsisanartifactofanincompletecorrelationtreatment. The comparison between the calculated and experimental spectrum confirmed D 2h struc- ture of the anion, as well as accuracy of the coupled-cluster and spin-flip geometries, frequenciesandnormalmodesoftheanionandthediradical. Densityfunctionalcalcu- lations(B3LYP)yieldedonlyaD 2h minimum,however,thequalityofthestructureand vibrational frequencies is poor, as follows from the comparison to the high-level wave function calculations and the calculated spectrum. The analysis of charge localization patterns and the performance of different functionals revealed that B3LYP underesti- matesthemagnitudeofvibronicinteractionsduetoself-interactionerror. x Chapter1 ImplementationoftheEffective FragmentPotentialmethod 1.1 Introduction Condensed phase encompasses the majority of matter on the Earth. Not surprisingly, manymolecularsystemsthataresignificantduetotheirbiological,electricalormechan- ical properties also exist in the condensed phase. Properties of individual molecules in suchsystemsareinfluencedbytheenvironment. Therefore,intermolecularinteractions playcrucialroleinchemicalandphysicalpropertiesofthecondensedphase. There are three distinct classes of methods for modeling the properties of the extended systems at molecular level: ab initio (QM), classical (MM) and hybrid (QM/MM)approaches. The ab initio methods treat a molecule of interest together with the neighboring molecules as one system at the same level of ab initio theory. The benefit of this approach is that a high accuracy can be achieved, with a proper choice of the ab ini- tiomethod. However,computationalcostsofmanyabinitiomethodsscalesteeplywith thesizeofthesystem[forinstance,asN 7 forCCSD(T) 1 ],whichmakesthefullabinitio treatmentoftheenvironmentaleffectsimpractical. New directions in the development of the ab initio methods with more favorable computationalcostsare: 1 • Linear scaling variants of ab initio methods 2 are based on local and/or density fitting(resolutionofidentity)approximations. Atpresent,theyareavailablefora fewabinitiomethods(HF 3 ,DFT 4 ,MP2 5 ,CC 6 ). • Separation techniques such as semi-classical divide and conquer 7 and frag- ment molecular orbitals 8 decompose the molecular system into fragments (e.g., molecules), which are described by individual and thus localized wavefunctions. Thewavefunctionsforeachfragmentarecomputedinthefieldofthefrozenwave- functionsofallotherfragmentsuntilself-consistencyisachieved. • Plane-wave DFT represents the wavefunction in terms of the plane waves 9 . It is suitable for systems with periodic boundary conditions (such as metals and crystals)andfeatureslinearscalingofthecomputationaleffort. On the other side of the theory spectrum, the whole molecular system can be mod- eledusingmolecularmechanics,thatisclassicallyusingempiricalforcefields 10–17 . The advantage of this approach is its efficiency - very large systems with many thousands of atoms can be modeled. However, molecular mechanics is incapable of modelingchemicalreactionsormolecularpropertiesdependingontheelectronicstruc- ture. Furthermore description of the intermolecular interactions is typically limited to parametrizedpairwisepotentialsandthuslacksaccuracy. In recent years, a number of hybrid methods combining the accuracy of ab initio andtheefficiencyofmolecularmechanicsbytreatingdifferentpartsofthesystemwith differentlevelsoftheoryhasbeendeveloped 18–21 . The hybrid methods can be subcategorized into continuum and discrete ones (Fig. 1.1). The continuum methods consider the environment to be uniform, as the namesuggests. Itseffectsonanindividualmoleculeisrepresentedwithapotentialthat depends only on the shape of the cavity the molecule allocates and the bulk properties 2 of the matter. The simplicity of continuum models makes them computationally effi- cient, but limit their accuracy. The discrete methods treat individual molecules of the environment explicitly, for instance by including all the neighbor molecules within a certain radius in the model. Unlike the fully ab initio methods, the molecules of the environmentaredescribedclassically, i.e., byamolecularmechanics(MM)forcefield. Thediscretemethodsaresuperiorinthesensethattheiraccuracycanbesystematically improvedby(1)increasingthenumberofneighboringmoleculesincludedinthemodel andby(2)usingahigherleveltheorytodescribethem. Overall, thecontinuummodels are preferred when accurate description of the bulk is essential. For example, a shift in the absorption spectrum of a molecule in a solvent can be well reproduced by a con- tinuum model. They break down when the environment is non-uniform (e.g., enzyme active center) or when relaxation is important due to the strong interaction between solvent and solute. Such cases call for discrete methods. However, to describe bulk behaviorthediscretemethodsneedtobecoupledwithasamplingtechnique(molecular dynamicsorMonteCarlo). Figure 1.1: Comparison of continuum (a) and discrete (b) models for a water molecule intheenvironmentofotherwatermolecules. 3 Depending on the way QM and MM parts are integrated with each other, the dis- creteQM/MMmethodscanbecategorizedasmechanicalembeddingandelectrostatics embedding 18 . MechanicalembeddingcomputestheMMparametersfromtheQMwave- function and then treats the whole system classically. Although mechanical embedding is simple and computationally cheap, the large drawback is that the perturbation of the QMwavefunctionduetoMMpartisignored. Electricalembeddingexplicitlyaccounts for the presence of the MM part by including the perturbative potential from MM into theHamiltonianoftheQMsystem. Thus,itprovidesbetterdescriptionoftheelectronic structureoftheQMpart,atthepriceoftheincreasedcomputationalcomplexity. One of the promising hybrid methods based on the electrostatic embedding is the effectivefragmentpotential(EFP)approach 22–33 . In EFP the molecular system is decomposed into an active region and a number of effective fragments. An example of such decomposition for uracil (active region) in water(effectivefragments)isdemonstratedinFig.1.2. Figure 1.2: Uracil in aqueous solution. Active region (uracil) is shown with a balls and sticksmodel. Effectivefragments(water)aredepictedbylines. Theactiveregionisapartofthemolecularsystemthatrequiresaccuratedescription, forinstance,becauseitparticipatesinachemicalreactionoritselectronicpropertiesare 4 ofinterest. Therefore,theactiveregionismodeledwithanabinitiomethod. Theeffec- tivefragmentsarethespectatormoleculesandthusdonotrequirehighaccuracyfortheir internal description. However, since the interaction between the active region and the effectivefragmentsaffectsthepropertiesoftheactiveregion,thisinteractionshouldalso be modeled with reasonable accuracy. An effective fragment in EFP is described semi- classically by a set of multipoles (with screening coefficients), polarizability points, dispersionpointsandexchange-repulsionparameters. The interaction between the effective fragments and the active region is accounted forbyaddingaone-electrontermtotheHamiltonianoftheactiveregion: H 0 =H +hp|V|qi, (1.1) whereV isthepotentialduetotheeffectivefragments. The expression for the inter-fragment interaction energy is derived by applying perturbation theory to the wavefunctions of two interacting molecules, as detailed in Ref.53. Withinthisapproach,thetotalenergyforfragmentsAandBis: E total =E 0 +E 1 +E 2 +... E 0 =E A +E B E 1 =h00|T|00i E 2 =− X mn h00|T|mnihmn|T|00i E 0 mn −E 0 00 , (1.2) whereE A andE B are energies of the individual isolated fragments, E n are the pertur- bativecorrections,T istheelectrostaticsoperator,E 0 mn istheenergyoftheexcitedstate mnforthetwonon-interactingfragments. 5 The first order perturbative correction E 1 is the electrostatic interaction. In EFP, it is described by a distributed multipole expansion (section 1.2). The second order perturbative correction consists of polarization and dispersion terms. The polarization component of EFP is accounted for by distributed polarizabilities (section 1.3). The dispersionisrepresentedsimilarlybyasetofdistributeddispersionpoints(section1.4). These terms are derived from the exact ab initio expressions and, therefore, their accu- racy can be systematically improved. The higher order terms in perturbative expansion are expected to have smaller values and are neglected in the current EFP formalism. One of the reasons of the EFP success 29 is that some of the higher order terms are believed to cancel out each other, thus omitting them does not affect the total inter- action energy substantially. For example, charge-transfer and exchange-induction are similar in magnitude, but have the opposite signs and thus cancel out each other. The exchange-dispersion partially cancels out with higher order dispersion terms (such as induceddipole-inducedquadrupole)thatarealsoomittedintheEFPformalism 29 . Perturbativetreatmentofintermolecularinteractionsdescribedaboveworkswellfor largedistances,butbreaksdownwhenmoleculesareclosetoeachotherandtheirelec- tronic densities overlap. To account for this effect, EFP includes a charge-penetration correction(section1.2)totheelectrostaticsandtheexchange-repulsioncomponent(sec- tion1.5)accountingforthePauliexclusionprinciple. Bothtermsarealsoderivedfrom the exact ab initio expressions. Thus, EFP intermolecular interaction consists of four components: electrostatics (with charge-penetration correction), polarization, disper- sionandexchange-repulsion. Theinteractionenergybetweentheeffectivefragmentsiscomputedas: E ef−ef =E elec +E pol +E disp +E ex−rep (1.3) 6 Since the active region is described by an ab initio method, the electrostatics and polarization interactions with the effective fragments are accounted for by additional one-electrontermsintheHamiltonian: H 0 =H + p V elec +V pol q (1.4) In the present EFP model, dispersion and exchange-repulsion interactions between active region and the effective fragments are treated similarly to the fragment-fragment interactionsasadditivecorrectionstothetotalenergy. OneoftheadvantagesofEFPisthatallofitstermshasbeenderivedfromtheexact ab initio expressions. Thus, the parameters necessary to describe an effective fragment can be obtained from the wavefunction of an isolated molecule. The parameters are computed only once for every fragment type and then can be reused for all consequent computations. This makes EFP method very efficient. Its computational cost isO(N 2 ) forthesystemwithN effectivefragments,plusthecostoftheabinitiocomputationfor theactiveregion,whichdoesnotdependonthenumberoffragments. AnotheradvantageofEFPisthat,unlikestandardQM/MMapproaches,itisapolar- ized force field. Therefore, the effective fragments not only polarize the active region, but also are polarized by it (polarized embedding), as well as by each other. Therefore, EFPdescribespolarizationeffectsmoreaccuratelythanstandardQM/MMapproaches. The third advantage is that the accuracy of EFP can be systematically improved by introducing higher order terms and cross terms in the perturbative expansion (e.g. exchange-dispersion), as well as higher order terms in the expansion of the individual components (e.g. hexapole-hexapole term for electrostatics, or induced dipole-induced quadrupoletermfordispersion). 7 The first implementation of EFP appeared in the GAMESS 34 quantum chemistry package, which is developed by the original authors of the EFP method 22 . By now EFP has been successfully applied for treating the environmental effects in a number of molecular systems 35–41 . In addition, due to the accurate description of the inter- fragmentinteractions,EFPhasbeensuccessfullyusedtostudytherelativestabilitiesof themolecularclustersandthedissociationenergyprofiles 27,28,42–46 . TheprimarygoalofthisworkwastoimplementtheEFPmethodinQ-Chem 47 . The benefitsofsuchimplementationare: (1)theextendedavailabilityoftheEFPmethod,(2) the ability to generate EFP fragments parameters with a higher level methods available in Q-Chem (such as CCSD) and (3) creating the base for future development of EFP in whichtheactiveregioncanbemodeledwithacorrelatedabinitiotheory a . The second goal was to apply the developed EFP implementation to the benzene clusters. Thesesystemsareinterestingduetotheirπ- π stackinginteractions,whichare similarinnaturetothoseinDNA 48 . The individual components of EFP are described in details in the four following sections. Section 1.6 describes the developed implementation of EFP in Q-Chem. Its applicationtobenzenetrimersandhigheroligomersispresentedinsection1.7. Thelast sectionofthechapterpresentssummaryoftheresultsandconclusion. a ThepresentGAMESSimplementationonlyenablesanon-correlateddescriptionoftheactiveregion (HForDFT). 8 1.2 Electrostatics Electrostatic component of the EFP energy accounts for Coulomb interactions. For molecular systems with hydrogen bonds or for strongly polar molecules, electrostatics hastheleadingcontributiontototalEFPenergy 49 . EFP electrostatics was devised to provide an accurate, yet compact and computa- tionally efficient description of the molecular electrostatic properties 22 . Buckingham 49 has shown that the Coulomb potential of a molecule can be expressed using multipole expansionatasinglespacialpointk: V elec k (x)= q k r kx + x,y,z X a μ k a F a (r kx )+ 1 3 x,y,z X a,b Θ k ab F ab (r kx )− 1 15 x,y,z X a,b,c Ω k abc F abc (r kx ), (1.5) whereq,μ,Θ,andΩarenetcharge,dipole,quadrupole,andoctopole,respectively,and F a , F ab , and F abc are the electric field, field gradient, and field second derivatives, due to the effective fragment. This expression converges to the exact electrostatic poten- tial when an infinite number of terms are included. However, it is known to converge slowly 50 and thus requires a large number of terms to achieve a reasonable accuracy, evenifsuchexpansionsareassignedtotheindividualatomsinthemolecule,asinMul- liken population analysis 51 . That makes a single point or atomic based representations ofthemolecularelectrostaticpotentialimpractical. In his work on distributed multipoles analysis 50,52–54 , Stone has shown that a better description of the Coulomb potential can be obtained by using many points as centers of the multipole expansions. In fact, for a system with a wavefunction expressed in the basisofN gaussians,theCoulombpotentialcanbedescribedexactlywithafinitenum- berofmultipoleslocatedatN 2 pointscorrespondingtothegaussianproductscenters 55 . However,thisrepresentationbecomesinefficientforsystemswithalargebasis. 9 Numericaltestsconductedonanumberofsystemshavedemonstratedthatanaccu- rate representation of the electrostatic potential can be achieved by using multipoles expansionaroundatomiccentersandbondmidpoints(thatis, thepointsknowntohave high electronic density) and truncating multipoles expansions after the octopoles 50,52 . ThisrepresentationhasbeenadoptedbyEFP 22,25 . Thus,electrostaticpotentialofaneffectivefragmentisspecifiedbycharges,dipoles, quadrupoles and octopoles placed at atomic centers and bond midpoints of a molecule. AnexampleofsuchexpansionforawatermoleculeisshowninFig.1.3. Figure1.3: Distributedmultipoleexpansionpointsforawatermolecule. The energy of electrostatic interaction for a system of effective fragments can be expressedas: E elec = 1 2 X A X B6=A X k∈A X l∈B E elec kl + X I∈A X J∈B E elec IJ + X I∈A X l∈B E elec Il ! , (1.6) where A and B are the fragments, k and l are the multipole expansion points, I and J arethenuclei. 10 The electrostatic interaction between two expansion points is computed using the followingequation 22 : E elec kl =E ch−ch kl +E ch−dip kl +E ch−quad kl + E ch−oct kl +E dip−dip kl +E dip−quad kl +E quad−quad kl (1.7) Thespecificequationsforcharge-charge(D.1),charge-dipole(D.2),charge-quadrupole (D.3), charge-octopole (D.4), dipole-dipole (D.5), dipole-quadrupole (D.6) and quadrupole-quadrupole(D.7)areprovidedinAppendixD. Nuclei-multipole and nuclei-nuclei terms are given by equations (1.8) and (1.9) respectively 22 . E elec Il =E nuc−ch Il +E nuc−dip Il +E nuc−quad Il +E nuc−oct Il (1.8) E elec IJ = Z I Z J R IJ (1.9) Equations for nuclei-charge (D.8), nuclei-dipole (D.9), nuclei-quadrupole (D.10) and nuclei-octopole(D.11)termscanbefoundinAppendixD. Electrostatic interaction between an effective fragment and the ab initio part is describedbyintroducingperturbationtoabinitioHamiltonian: H =H 0 +V (1.10) Theperturbationisaone-electrontermcomputedasasumofcontributionsfromallthe expansionpointsandnucleiofalltheeffectivefragments: V = X p X q X A X k∈A d pq p V elec k q + X I∈A p Z I R q ! , (1.11) 11 where A is an effective fragment, I is a nucleus, p and q are atomic orbitals of the ab initiopart,Z I isnucleuscharge,d pq isanelementoftheatomicdensitymatrixandV elec k iselectrostaticpotentialfromanexpansionpoint. A one-electron contribution to the Hamiltonian from an individual expansion point consists of 4 terms, each of which originates from the electrostatic potential of the cor- respondingmultipole 22 : p|V elec k |q = p q k R q + p P x,y,z a μ k a a R 3 q + * p P x,y,z a,b Θ k ab (3ab−R 2 δ ab ) 3R 5 q + + * p P x,y,z a,b,c Ω k abc (5abc−R 2 (aδ bc +bδ ac +cδ ab )) 5R 7 q + , (1.12) whereR,a,b andc are the distance and distance components between the electron and theexpansionpointk. Electrostaticsparameters The parametersq,μ,Θ,Ω of the effective fragment electrostatics need to be computed only once for every fragment type. These multipoles are obtained from an ab initio electronicdensityoftheindividualmoleculeusingDMA(distributedmultipolesanaly- sis) procedure described by Stone 50,52 . First, a molecular wavefunction expressed in a basis ofN gaussians is computed with an ab initio method. Then spherical multipoles expansions are calculated at the centers of the gaussian basis function products. These expansions are finite and include multipoles up to L 1 +L 2 order where L 1 and L 2 are angularmomentsofthe two gaussians. Atthenext step, thesemultipolesare translated totheEFPexpansionpoints(atomiccentersandbondmidpoints)asdescribedinRef.50 12 and truncated after the4 th term. Thus obtained set of spherical multipoles is then con- vertedtotracelessBuckinghammultipolesusingtheequationsprovidedinAppendixE. Chargepenetration ModelofEFPelectrostaticsdescribedaboveworkswellwhenelectronicdensitiesofthe effective fragments do not overlap. However, for many molecular systems the overlap issubstantialandisknownascharge-penetrationeffect. Asketchdemonstratingcharge penetrationfortwowatermoleculesisshowninFig.1.4 Figure1.4: Sketchofthechargepenetrationphenomenonfortwowatermolecules. The magnitude of this effect is commonly around 15% of the total electrostatic energyandforsomesystemsisaslargeas200% 27 . ItwasdemonstratedinRef.24that charge penetration can be efficiently accounted for by introducing a screening function asaprefactortotheCoulombpotential: V elec μk (x)→ 1−e −α μk r x−μk V elec μk (1.13) 13 This screening is applied only to the potential originating from charges representing electronic density, and not to the potential of higher multipoles, which represents elec- tronicdensityvariations 27 . Thereforeonlythecharge-chargetermismodifiedtoaccount forthecharge-penetration 24 : E ch−ch kl = 1−(1+ α k R kl 2 )e −α k R kl q k q l r kl ifα k =α l 1− α 2 l α 2 l −α 2 k e −α k R kl − α 2 k α 2 k −α 2 l e −α l R kl q k q l r kl ifα k 6=α l (1.14) Electrostaticsscreeningparameters Parameters a and b of the screening function are computed using a fitting procedure described in Ref. 24 once for every type of the fragment. This procedure optimizes the screening parameters so that the sum of differences between the ab initio electrostatic potentialandtheEFPelectrostaticpotentialwithscreeningisminimized: Δ= X p V Elec abinitio (p)−V Elec EFP (p) , (1.15) wherethesumistakenoverallthepointspofthegrid. Thegridisconstructedsuchthat itincludesthepointswithsubstantialelectronicdensityandexcludesthepointscloseto nuclei and bond midpoints, which are inaccessible to multipoles of another fragments. AprojectionofsuchgridforawatermoleculeisshowninFig.1.5. 1.3 Polarization Polarizationaccountsfortheintramolecularchargeredistributionundertheinfluenceof an external electric field. It is the major component of many-body interactions, which 14 Figure1.5: Theshapeofthegrid(monotonousgrayarea)forfittingthescreeningparam- etersofawatermolecule. are important for proper description of cooperative molecular behavior. Since polariza- tioncannotbeparametrizedbytwo-bodypotentials,itisoftendisregardedinmolecular mechanics models. However, for hydrogen bonded systems, up to 20% of the total intermolecularinteractionenergycomesfrompolarization 22 . Theexactexpressionforthepolarizationenergyappearsasthesecond-ordertermof along-rangeperturbationtheory 53 : E pol =− X n6=0 h00|V|0nih0n|V|00i E n −E 0 , (1.16) wherenistheexcitedelectronicstateandV istheperturbingelectrostaticpotential. ExpandingelectrostaticpotentialV in 1 R n seriesinequation(1.16)andthencom- biningtermswiththesameelectricfieldderivativesyieldsthefollowingequation: E pol =− 1 2 α ab F a F b − 1 6 β abc F a F b F c − 1 3 A a,bc F a F bc − 1 6 B ab,dc F a F b F cd − 1 6 C ab,dc F ab F cd +..., (1.17) 15 whereF a andF ab areelectricfieldderivatives,α isthedipolepolarizabilitytensor,β is thedipolehyperpolarizability,A,B andC arequadrupolepolarizabilitytensors. Inordertoachieveabetterconvergence,EFPusesdistributedpolarizabilitiesplaced at the centers of valence LMOs. This allows one to truncate the expansion (1.17) after the first term and still retain a reasonable accuracy 22 . It should be noted that unlike the total molecular polarizability tensor, which is symmetric, the distributed polarizability tensors are asymmetric and thus all of 9 tensor components are required. Polarizability pointsforawatermoleculeareshowninFig.1.6. Figure1.6: Distributedpolarizabilitiesforawatermolecule. Letusdiscusshowpolarizationenergyiscomputedforanindividualeffectivefrag- ment in an external force field. First, induced dipoles for every polarizability point are computed: μ k =α k F k , (1.18) where μ k , α k and F k are the induced dipole, the polarizability tensor and the external forcefieldatpointk. 16 Theninduceddipolesareusedtocomputethepolarizationenergy: E pol = X k μ k F k (1.19) Now, let us consider a system with many effective fragments and an ab initio part. In this case, the electric field acting on an effective fragment depends on: multipoles andnucleiofothereffectivefragments,induceddipolesofothereffectivefragmentsand electronicdensityandnucleioftheabinitiopart: F k = X B6=A X i∈B F mult i (k)+ X I∈B F nuc I (k)+ X j∈B F ind j (k) ! +F abinitio elec (k)+F abinitio nuc (k), (1.20) whereF k isthetotalelectricfieldatpolarizabilitypointk oftheeffectivefragmentA,j isapolarizabilitypoint,iisamultipoleexpansionpoint,I isanucleus,F mult i (k)isthe field due to the multipolei,F nuc I (k) is the field from the nucleusI,F ind j (k) is the field due to the induced dipole j, F abinitio elec (k) and F abinitio nuc (k) are the forces from the electrondensityandnucleioftheabinitiopart. The dependence of the total electric field on other induced dipoles and ab initio density gives rise to two complications in computing the polarization energy discussed below. The first complication arises from the fact that an induced dipole depends on the values of all other induced dipoles in all other fragments in according with equations (1.18)and(1.20). Therefore,theymustbecomputedself-consistently. 17 Anothercomplicationoriginatesfromthedependenceoftheinduceddipolesonthe abinitioelectrondensity,which,inturn,isaffectedbythefieldcreatedbytheseinduced dipolesthroughaoneelectroncontributiontotheHamiltonian: H ind = X A X k∈A P x,y,z a (μ a k + ¯ μ a k )a R 3 , (1.21) whereR,a,b andc are the distance and distance components between the electron and the polarizability point k, μ a k and ¯ μ a k are components of the induced dipole and conju- gated induced dipole, respectively. Thus, when the ab initio part is present, conjugated induceddipolesshouldalsobecomputed: ¯ μ k =α T k F k , (1.22) whereα T k andF k are the transposed polarizability tensor and the external force field at the polarization point k. In the following discussion it is implied that the conjugated induceddipolesarealwayscomputedtogetherwiththeregularinduceddipoles. As a result of these self-consistency requirements, the total polarization energy is computedwithatwoleveliterativeprocedureillustratedinFig.1.7. Theobjectiveofthe higherlevelistoconvergethewavefunction. Thelowerlevelistaskedwithconverging theinduceddipolesforagivenfixedwavefunction. Inthebeginning,theinduceddipolesareinitializedto0andtheelectronicdensityis obtainedthroughaguessprocedure(suchasageneralizedWolfsberg-Helmholtzproce- dure 56 ). Ateveryiterationofthehigherlevel,thewavefunctionisupdatedbasedonthevalues oftheinduceddipolesonthepreviousstep. Thentheforcefromthisnewwavefunction acting on the polarizability points is recomputed. The convergence criteria is that the 18 Is μ converged? μ = 0 Ψ = guess() Compute Ψ Compute F Ψ Compute F μ Compute μ No Is Ψ converged? Yes Exit No Yes Figure1.7: Computingtheinduceddipolesthroughthetwolevelself-consistentiterative procedure. HereΨistheabinitiowavefunction,μisasetofinduceddipoles,F Ψ isthe forceduetotheabinitiowavefunctionandF μ istheforceduetotheinduceddipoles. wavefunction parameters (molecular orbital coefficients) are within a predefined range fromtheirvaluesonthepreviousiteration. Atthelowerlevel,thevaluesoftheinduceddipolesarecomputedateveryiteration based on the electric force from the ab initio part (which does not change during the lower-level iterations), as well as the forces from multipoles, nuclei and the force from theinduceddipolesofthepreviouslower-leveliteration. Theconvergencecriteriaisthat thedifferencebetweenthenewinduceddipolesandthosefromthepreviousiterationis withinapredefinedthreshold. Thus, when the lower-level iterative procedure exits the induced dipoles are self- consistent and are consistent with a fixed ab initio wavefunction. Convergence of the 19 higher-level iterative procedure results in the induced dipoles and the ab initio wave- functionthatareself-consistentandconsistentwitheachother. Ifthemolecularsystem doesnotcontaintheabinitiopart,thenonlythelower-levelprocedureisrequired. Once the values of the induced dipoles are obtained, the total energy of the system canbecomputedas 22 : E =E other − 1 2 X A X k∈A μ k (F mult k +F nuc k )+ 1 2 X A X k∈A (μ k + ¯ μ k )F abinitio elec k , (1.23) where E other is the sum of energies of the ab initio part and other components of EFP, Aisaneffectivefragment,k isapolarizabilitypoint,μ k and ¯ μ k areinduceddipolesand conjugatedinduceddipoles,F mult k istheforceduetothemultipoleexpansionpointsof otherfragments,F nuc k istheforceduetotheallothernuclei,F abinitio elec k istheforcedue totheabinitioelectronicdensity(F abinitio elec k = D Ψ ˆ f elec k Ψ E ). Inequation(1.23),thesecondtermisthetotalpolarizationenergyofalltheeffective fragments. However,thesecondandthethirdtermdonotadduptothetotalpolarization energy because a part of the later is implicitly included in the ab initio energy (due to theeffectoftheinduceddipolesonthewavefunction). Computationalcostofthepolarizationcomponentconsistsofthecostofcomputing theinduceddipolesandtheenergy. ThelatterisatrivialstepwiththecostofO(N 2 )for N effectivefragments. Thus,theformerdefinestheoverallcomputationalcost. At the higher level, the number of iterations h is similar to that for the ab initio systemwithoutEFP,whileatthelowerleveltheconvergenceisachievedwithinbounded number of iterations. Every lower-level iteration scales as O(N 2 ), thus, the overall computationalcostisO(hN 2 ). If system does not have the ab initio part, then only the lower-level iterative proce- dureisrequiredandcomputationalcostisO(N 2 ). 20 Polarizationparameters The polarization parameters in the EFP model are LMOs and polarizability tensors. Both can be obtained from an ab initio calculation once for every type of the fragment andthenreusedinallconsequentcalculations. LMOsareobtainedwiththeBoyslocal- izationprocedurebasedontheabinitioelectronicdensity 57,58 . Polarizabilities are obtained based on the definition (1.24) as a derivative of the dipolemoment(μ)w.r.t. electricfield(F)usingthefinitedifferencesmethod. α ab = ∂μ a ∂F b , (1.24) Inthisapproach,adipolemomentofanLMOisfirstcomputedwithoutanelectricfield. Then,bymeansof3additionalabinitiocomputations,theLMOanditsdipolemoment arecomputedinthepresenceofasmallelectricfieldactingalongoneofthecoordinate axis. BasedonLMOdipolemoments,thepolarizabilitytensorcanbecomputedas: α ab = μ a (F b )−μ a (0) F b (1.25) 1.4 Dispersion For many molecular systems dispersion is the main attractive component of the inter- molecularinteractionenergy. Anaccuratedescriptionofdispersioninteractionisadiffi- culttaskevenforabinitiomethods—inclusionoftheelectroncorrelationathighlevel is required. An additional challenge for EFP is that its dispersion description must be computationallyefficient. 21 The exact expression (1.26) for the dispersion energy appears in the second order termofalong-rangeperturbationtheory. E disp =− X n6=0 X m6=0 h00|V|mnihmn|V|00i E A m +E B n −E A 0 −E B 0 , (1.26) whereAandB aretwointeractionmolecules,mandnaretheirexcitedelectronicstates andV istheperturbingelectrostaticpotential. By expanding electrostatic potential V in 1 R n series dispersion energy can be expressedinthesocalledLondonseries: E disp = C 6 R 6 + C 8 R 8 + C 10 R 10 +..., (1.27) where theC n coefficients correspond to induced dipole—induced dipole (C 6 ), induced dipole—induced quadrupole (C 8 ), induced quadrupole—induced quadrupole and induced dipole—induced octopole (C 10 ) etc. interactions, Thus, dispersion energy can bealternativelyexpressedasasum: E disp =E ind dip−ind dip +E ind dip−ind quad +E ind quad−ind quad +E ind dip−ind oct +... (1.28) A better convergence can be obtained by using a distributed dispersion model with expansion points located at localized molecular orbitals centroids 59 . Expansion points forawatermoleculearedemonstratedinFig.1.8. When using distributed dispersion, a reasonable accuracy can be achieved by trun- cating the expansion after the first term and approximating the rest as 1/3 of the this 22 Figure1.8: Distributeddispersionpointsforawatermolecule. term 59 . Induced dipole-induced dipole interaction energy between two induced dipoles kandlcanbecomputedbyusingthefollowingequation: E ind dip−ind dip kl = x,y,z X abcd T kl ab T kl cd Z ∞ 0 α k ac (iν)α l bd (iν)dν, (1.29) whereα k ac andα l bd aredynamicpolarizabilitytensors,iν isthetheimaginaryfrequency, T kl ab andT kl cd aretheelectrostatictensorsofthesecondrank,i.e.,T kl ab = 3ab−R 2 kl δ ab R 5 kl . Ingeneral,theinteractionenergydescribedby(1.29)isanisotropic,distancedepen- dent and computationally expensive 26 . The approximation used by EFP includes only thetraceofpolarizabilitytensor,whichisanisotropicanddistanceindependent. Thetotaldispersionenergyforasystemofeffectivefragmentsis: E disp = 4 3 X A X B X k∈A X l∈B C kl 6 R 6 kl , (1.30) 23 whereAandB areeffectivefragments,k andlareLMOs,C kl 6 isintermoleculardisper- sioncoefficient,andR 6 kl isthedistancebetweentwoLMOcentroids. Intermolecular dispersion coefficients are computed using the 12 point Gauss- Legendreintegrationformula: C kl 6 = 12 X i=1 w i 2ν 0 (1−t i ) 2 ¯ α k (iν i )¯ α l (iν i ), (1.31) where ¯ α k (iν i ) and ¯ α l (iν i ) are1/3 of the traces of the corresponding polarizability ten- sors. Dispersiondamping For small distances between effective fragments dispersion coefficients must be cor- rectedtotakeintoaccountthechargepenetrationeffect. IntheoriginalEFPmodel,the Tang—Toenniesformula 60 withparameterb=1.5isused: C kl 6 → 1−e −bR ∞ X k=0 (bR) k k! ! C kl 6 (1.32) However,recentstudiesrevealedthatthisformulaover-dampsthedispersionandamore appropriatevalueforbliesbetween2and3andissystemdependent 29 . Asanalternative, Ref.29suggestedtouseadampingfunctiondependentontheelectronicdensityoverlap S kl betweenLMOsk andl: C kl 6 → 1−|S kl | 2 (1−2ln|S kl |+2ln 2 |S kl |) C kl 6 (1.33) 24 Dispersionparameters As follows from equations (1.30) and (1.31), the dispersion parameters for an effective fragment are the coordinates of LMO centroids and sets of dynamic polarizability ten- sor traces for 12 frequencies for each LMO. These parameters can be computed once for every type of an effective fragment and then reused in EFP calculations. Dynamic polarizabilities can be computed using dynamic time-dependent Hartree-Fock (TDHF) ordynamictime-dependentdensityfunctionaltheory(TDDFT)asdescribedinRef.26. 1.5 Exchange-repulsion Exchange-repulsion originates from the Pauli exclusion principle, which states that for two identical fermions the wavefunction must be anti-symmetric. In ab initio methods suchasHartree-FockandDFT,itisaccountedforbytheexchangeoperator. Inclassical force fields, exchange-repulsion is introduced as a positive (repulsive) term, e.g., 1 r 12 in the Lennard-Jones potential. However, since the Pauli exclusion principle is intrin- sically quantum, its classical descriptions are not accurate. EFP uses a wavefunction- based formalism to account for the electron exchange. Thus exchange-repulsion is the onlynon-classicalcomponentofEFPandtheonlyonethatisrepulsive. The EFP exchange-repulsion energy term has been derived from the exact equation fortheexchange-repulsionenergy 23,61 : E exch = D Ψ A Ψ B | ˆ AH AB |Ψ A Ψ B | E D Ψ A Ψ B | ˆ AΨ A Ψ B E −hΨ A Ψ A |V AB |Ψ B Ψ B i−E A −E B , (1.34) 25 whereA andB are two molecules described byΨ A andΨ B wavefunctions,V AB is the intermolecularCoulombinteraction,and ˆ Aistheantisymmetrizationoperator: ˆ A=P 0 −P 1 +P 2 −P 3 +..., (1.35) where P i is a permutation of i electron pairs consisting of one electron from each molecule. Truncating sequence (1.35) after the second term, applying an infinite basis set approximation and a spherical gaussian overlap approximation as shown in Ref. 61 leadstotheexpressionfortheexchange-repulsionenergyintermsoflocalizedmolecu- larorbitals: E exch = 1 2 X A X B6=A X i∈A X j∈B E exch ij (1.36) E exch ij =−4 r −2ln|S ij | π S 2 ij R ij −2S ij X k∈A F A ik S kj + X l∈B F B jl S il −2T ij ! +2S 2 ij X J∈B −Z J R −1 iJ +2 X l∈B R −1 il X I∈A −Z I R −1 Ij +2 X k∈A R −1 kj −R −1 ij ! , (1.37) whereA andB are the effective fragments,i,j,k andl are LMOs,I andJ are nuclei, S and T are intermolecular overlap and kinetic energy integrals, F is the Fock matrix element. ThepositionsofLMOcentroidsforawatermoleculeareshowninFig.1.9. Formula for theE exch ij involves intermolecular overlap and kinetic energy integrals, whichareexpensivetocompute. Inaddition, sinceequation(1.37)isderivedwithinan 26 Figure1.9: Theexchange-repulsionLMOcentroidsforawatermolecule. infinite basis set approximation, it requires a reasonably large basis set to be accurate. These factors make exchange-repulsion the most computationally expensive part of the EFPenergycalculationsofmoderatelysizedsystems. Large systems require additional considerations. Equation (1.36) contains a sum over all the fragment pairs and its computational cost formally scales as O(N 2 ) with the number of effective fragments N. However, exchange-repulsion is a short range interaction —- the overlap and kinetic energy integrals of LMOs vanish very fast with the interfragment distance. Therefore, by employing a distance based screening, the number of overlap and kinetic energy integrals scales as O(N). As a result, for large systems exchange-repulsion becomes less computationally expensive than long range componentsofEFP(suchasCoulombinteractions). Exchange-repulsionparameters Theparametersrequiredforequation(1.37)aretheFockmatrixelementsandtheLMOs (orbitalcoefficients,basisandcentroids). Theycanbecomputedonceforeverytypeof 27 thefragmentbyperformingaHartree-Fockcalculationonthemoleculefollowedbyan anorbitallocalizationprocedure(suchasBoysLMO) 23 . 1.6 ImplementationofEFPinQ-Chem ThissectiondescribesthenewimplementationoftheEFPmethodintheQ-Chemelec- tronicstructurepackage. As detailed in section 1.1, within EFP formalism the molecular system is decom- posedintoanactiveregionandanumberofeffectivefragments. Sincetheactiveregion istreatedbyanabinitiomethod,EFPrequiresaninterfacewithageneral-purposequan- tum chemistry package. Our implementation was developed using C++ 62 as a module withinthestateof-the-artquantumchemistrypackageQ-Chem 47 . Thestandardtemplate library(STL)wasextensivelyusedthroughouttheimplementation. The developed code consists of the two independent parts: EFP energy and EFP parameters. The former is responsible for computing energy of the molecular system, whereasthelatterallowsonetoobtaineffectivefragmentparameters. ComputingEFPenergy The EFP energy part of the implementation allows one to compute total energy of a system comprised of an optional active region and a number of effective fragments. The total energy includes contributions from all the EFP components described in the fourprevioussections: electrostatics(withcharge-penetrationcorrection),polarization, dispersion,exchange-repulsion. Asimplifiedhigh-levelclassdiagramfortheEFPenergycodeisshowninFig.1.10 (a dashed arrow represents a dependence). The main class EFP is responsible for man- aging the interactions between other classes and providing the interface to the rest of 28 the Q-Chem package. Classes Electrostatics, Polarization, Dispersion and Exchange- repulsionencapsulatethefunctionalityrequiredforthecorrespondingindividualenergy components. Class Polarization has the dependence on the class Electrostatics because induced dipoles depend on the field created by multipoles. The responsibility of the ConfigReader class is reading the effective fragment parameters and configuration of themolecularsystem,andstoringthemintheinternalrepresentationusingclassesFrag- mentParamsandFragment. D i s p e r s i o n E F P E l e c t r o s t a t i c s P o l a r i z a t i o n E x c h a n g e R e p u l s i o n C o n f i g R e a d e r F r a g m e n t P a r a m s F r a g m e n t Figure1.10: ClassdiagramfortheEFPenergypartoftheimplementation. FragmentParams std::vector< double > er_fock_matrix std::vector< PolarizationPoint > polarizabilities std::vector< XYZPoint > er_lmos std::vector< std::vector< double > > er_wavefunction std::vector< DispersionPoint > dispersions std::string er_basis std::vector< EFPAtom > atoms std::vector< MultipolePoint > multipoles Figure1.11: TheinternalstructureoftheFragmentParamsclass. 29 The internal structure for FragmentParams class is shown in Fig. 1.11. The labels above the arrows indicate the name of the data member. This class encapsulates all the parameters for a given EFP fragment type: atoms, multipoles, dispersion points, polarizabilitiesandparametersforexchange-repulsion(whichstartwither prefix). TheclassFragmentrepresentsanindividualeffectivefragmentinthemolecularsys- tem. It stores the fragment type, its position and orientation. The position is specified in terms of translation in Cartesian coordinates (x, y, z) relative to the origin of the ab initio system. The orientation is specified by the Euler angles (α, β, γ) in terms of the 3 consequent rotations, as explained in Fig. 1.12. The rotation matrix corresponding to theseanglesisgiveninFig.1.13. N x y z Z X Y α β γ Figure1.12: RepresentingtheorientationwiththeEulerangles. xyzisthefixedsystem, XYZ is the rotated system, N is the intersection between the xy and XY planes called the line of nodes. Transformation between the xyz and XYZ systems is given by α rotation in the xy plane, followed by β rotation in the new y’z’ plane (around the N axis),followedbyγ rotationinthenewx”y”plane. 30 R= cosαcosγ−sinαcosβ sinγ −cosαsinγ−sinαcosβ cosγ sinβ sinα sinαcosγ+cosαcosβ sinγ −sinαsinγ+cosαcosβ cosγ −sinβ cosα sinβ sinγ sinβ cosγ cosβ Figure1.13: TherotationmatrixcorrespondingtotheEuleranglesα,β,γ. Thepositionx 0 i ,y 0 i ,z 0 i ofanatomiofaneffectivefragmentcanbecomputedas: x 0 i y 0 i z 0 i = R x i y i z i + x y z , (1.38) wherex i ,y i ,z i are the coordinates of the atom in the fragment frame, R is the rotation matrixdefinedinFig.1.13,andx,y,z arethefragmentcoordinates. TheflowofcontrolbetweenQ-ChemandtheEFPmoduleispresentedinFig.1.14. The ReadParams and ReadFragments calls instruct the EFP module to read the defini- tions of the fragment types and the molecular system configuration from the Q-Chem standard input (discussed in the next subsection). The Init call checks the effective fragments parameters and the configuration for consistency, creates and initializes the instances of Electrostatics, Polarization, Dispersion and Exchange-Repulsion classes. TheUpdateRepresentationmethodiscalledeverytimewhenthegeometryofthemolec- ularsystemchanges. Itperformsrotationandtranslationofatoms,multipoles,polariza- tion tensors, dispersion points and LMOs in accordance with the coordinates and ori- entation for each individual fragment. The NuclearEnergy function computes the EFP energy components, which are independent from the wavefunction of ab initio system (inter-fragmentelectrostatics,dispersionandexchange-repulsion). 31 Iftheactiveregionispresent,thenforeveryiterationoftheSCFprocedureQ-Chem callstheHamiltonianmethodandtheWFDependentEnergyfunction. TheHamiltonian method updates the one-electron term of the active region Hamiltonian in accrodance with the equations (1.10) and (1.21). After updating the wavefunction accordingly, Q- Chem calls WFDependentEnergy to compute the wavefunction-dependent part of the EFP energy (polarization). WFDependentEnergy takes the electronic density of the active region as a argument. If the active region is absent, Q-Chem calls WFDepen- dentEnergyonlyoncewithazeroelectronicdensityasanargument. Finally, the PrintSummary method is called to print the individual energy compo- nentsfortheinter-fragmentEFPenergy. Q-Chem EFP ReadParams() ReadFragments() Init() UpdateRepresentation() NuclearEnergy() Hamiltonian() WFDependentEnergy() energy energy for every SCF cycle PrintSummary() hamiltonian for every new geometry Figure1.14: TheflowofcontrolbetweenQ-ChemandEFPenergymodule. 32 The correctness of the EFP energy implementation in Q-Chem was verified by performing the energy computation on a few dozens of examples and comparing the resultswiththeGAMESSimplementationdevelopedbytheoriginalauthorsoftheEFP method. The examples used for testing include systems composed of water, ammonia andbenzenemoleculesindifferentconfigurations. EFPinputformat TheQ-Chemstandardinputconsistsofanumberofsectionsdelimitedby $keywordand $end. Themandatorysectionsare$molecule,whichdescribesthegeometryofabinitio part, and $rem section, which specifies parameters of the ab initio computation. Other sectionsarespecifictodifferentQ-Chemmodules. The EFP code extends the Q-Chem input with 2 new sections: $efp params and $efp fragments, which describe the parameters of the EFP fragment types and define themolecularsystembyprovidingfragments’positionsandorientations,respectively. In the EFP input sections, all coordinates are in Angstroms, unless Q-Chem’s INPUT BOHR keyword is set to TRUE, in which case all the coordinates are in Bohr, all angles are in radians and other values are in atomic units. All keywords used in the EFPinputsectionarecaseinsensitive. Fig. 1.15 presentswhenhe format for the $efp params section. It consists of the parameters’definitionsforoneormorefragmenttypes. Adefinitionstartsbyproviding the name of the fragment type, listing all atoms and their coordinates. It is followed bythelistofallthemultipoleexpansionpoints(mult). Eachmultipoleexpansionpoint is defined by its coordinate, multipoles and a damping parameter (DAMP). Multipoles withdifferentangularmomentumarespecifiedontheseparatelinesinanarbitraryorder. In thier definition C is charge, D , Q and O are dipole, quadrupole and octupole components,respectively. Allmultipolesareoptional,butwhenprovided,thenalltheir 33 components should be specified in the order demonstrated in Fig. 1.15. The damping parameterisoptionalandshouldbepreceededbythecdampkeyword. Next in the fragment definition are the polarization points (pol). They are specified bythecoordinatesandcomponentsofthepolarizationtensorP . Thedispersionpoints(disp)aredescribedsimilarlybythecoordinatesandthetraces ofthedynamicpolarizabilitytensorsC at12differentfrequencies. The last subsection describes the parameters required for the exchange-repulsion: er basisspecifiesthebasisused,er fockprovidestheFockmatrix,er lmosspecifiesthe LMO centroids and er wavefunction defines each LMO in terms of the atomic orbitals. The subsections for different EFP components are optional and can appear in an arbi- traryorder. Fig. 1.16 presents the format for the $efp fragments section. Every effective frag- mentinthemolecularsystemisdefinedbythefragmenttype,followedbyitsCartesian coordinates (X, Y, Z) specifying the translation relative to the ab initio system origin and Euler angles specifying the fragment orientation (α, β, γ). If the orientation is not provided, the default values of the angles are set to 0. The fragment type should match oneofthefragmentsinthe$efp paramssection. The active region is specified in the$molecule section the same way as an ab initio system is specified in Q-Chem 47 . In order to define a molecular system composed of the effective fragments only, the keyword EFP FRAGMENTS ONLY should be set to TRUE inthe$remsectionoftheQ-Chemstandardinput b . A complete example of the Q-Chem standard input for the system composed of 3 watermolecules(oneasanactiveregionandtwoothersaseffectivefragments)isgiven inAppendixF. b Since the $molecule section is mandatory in Q-Chem, even if the active region is disabled it should still contain at least one atom, for instance, He. In the nutshell, assigning EFP FRAGMENTS ONLY keywordtoTRUE simplydisablesanyinteractionbetweentheabinitiopartandtheeffectivefragments. 34 $efp params fragment <fragment type 1> ATOM X Y Z ... mult X Y Z [cdamp DAMP] [C] [DX DY DZ] [QXX QXY QXZ QYY QYZ QZZ] [OXXX OXXY OXXZ OXYY OXYZ OXZZ OYYY OYYZ OYZZ OZZZ] ... pol X Y Z PXX PXY PXZ PYX PYY PYZ PZX PZY PZZ ... disp X Y Z S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 ... er basis BASIS er fock F1 F2 F3 F4 ... er lmo X Y Z ... er wavefunction 1 C1 C2 C3 C4 ... 2 C1 C2 C3 C4 ... ... fragment <fragment type 2> ... $end Figure1.15: Theformatforthe$efp paramssection. $efp fragments <fragment type> X Y Z [α β γ] ... $end Figure1.16: Theformatforthe$efp fragmentssection. ComputingEFPparameters TheEFPparameterspartofthecodecomputestheelectrostatics,charge-penetrationand polarization parameters for an isolated molecule. The parameters need to be computed 35 onlyonceforeveryfragmenttypeandthencanbereusedinallsubsequentEFPenergy computations. Since the EFP model has been derived from the ab initio theory, the first step in computing the EFP parameters for a molecule is to obtain its electronic density. This can be achieved by means of any ab initio method available in Q-Chem that supports explicitelectronicdensity(forinstance: HF,DFT,CCSD). From the electronic density, the electrostatic parameters (multipoles) are computed usingthedistributedmultipoleanalysisprocedureasdescribedinsection1.2. Themul- tipole expansion is computed up to octupoles. The points for the expansion are the nucleiandthebondmidpoints. The screening parameters accounting for the charge-penetration are computed by the fitting procedure, which minimizes the difference between the exact electrostatic potentialandpotentialfromthemultipolesexpansionpointsmultipliedbythescreening function, as detailed in the section 1.2. The optimization of screening parameters is basedontheNewton-Raphsonalgorithm 63 . The polarization parameters are computed using the finite difference method based on 4 ab initio computations: one with no external field and 3 with a small (by default 10 −5 V/m) electrostatic field along one of the coordinate axis, as described in section 1.3. In each computation, the localized molecular orbitals are obtained with the Boys localizationprocedure 57,58 . The computation of the dispersion and exchange-repulsion parameters have not yet beenimplementedinQ-Chem. Theycanbeobtainedfromthequantumchemistrypack- age GAMESS 34 , using the code developed by the original authors of the EFP method. Usingitsmakefpjobtype,alltheEFPparameterscanbeproducedintheGAMESSfor- mat. AnexampleoftheGAMESSinputfilegeneratingtheEFPparametersforbenzene and the output produced are shown in Appendix A. To facilitate the conversion from 36 the GAMESS format to the Q-Chem format, a converter script (Appendix C) has been developedinPython 64 . IttakesafilewiththeGAMESSEFPparametersasacommand line argument and prints out the parameters in the Q-Chem format. The conversions performedinternallybythescriptarethesyntaxconversion,unitconversion(fromBohr to ˚ A) and redundancy elimination (replacing a dynamic polarizability tensor with its trace). The Q-Chem keywords affecting the EFP parameters computation are provided in AppendixG. 1.7 EFPapplication: Beyondbenzenedimer Aromatic π − π interactions have the major influence on the structure and properties of many biologically important molecular systems 65 . They govern nucleotide stacking patternsinDNAandRNA 48 ,controltheternarystructuresofproteins 66 ,playanimpor- tant role in the intercalation of drugs with DNA 67 , affect the properties of the techno- logically important polymers aramids 68 , impact the crystal structure of many aromatic compounds 69 ,andstabilizehost—guestcomplexes 70,71 . Studyingtheπ−πinteractionsisamajorchallengeforbothexperimentandtheory 72 . Theexperimentalinvestigationsaredifficultduetothecomplexenvironmentsinwhich theseinteractionsoccur(proteins,DNA).Theweaknessofindividualπ−π interactions and the flatness of the PES along the coordinates defining the relative position of the aromatic rings make it difficult to compute even for small model systems in the gas phase 73–75 . Thetheoreticalstudiesofπ−π interactions,asofanynon-covalentinteractions,are challengingduetotheweaknessoftheseinteractionsrelativetotypicalcovalentinterac- tions. Their proper treatment requires high accuracy ab initio methods [e.g. CCSD(T)] 37 and large basis sets with diffuse functions [e.g. aug-cc-pVTZ] 76 , which both contribute tothehighcomputationalcomplexity. Atpresent,highaccuracyabinitiocomputations arefeasibleonlyforsystemswithafewπ−π interactingaromaticrings. However,innatureseveralinteractingπ−π moietiesareusuallypresent,e.g.,many nucleotidesareresponsibleforthesecondarysecondarystructureofDNA 48 .Therefore, itisimportanttounderstandwhethermultipleπ−π interactionscanbetreatedasasum ofthepair-wiseinteractionsbetweenneighboringaromaticringsortherearesignificant non-additivemany-bodyinteractioncomponents. Theoreticalstudiesofmultipleπ−π interactionshavebeensofarlimitedduetothe computational complexity reasons discussed above 77–80 . An interesting approach was taken in the work by Tauer and Sherrill 78 who studied benzene trimers and tetramers usingMP2withaspeciallyselectedbasisset. MP2isknowntooverestimatethebinding energies. For instance, MP2 with a large basis set (aug-cc-pVQZ) overestimates the dissociation energy of the benezene dimer by 70% 76 . In contrast, small basis sets are knowntounderestimatebindingenergies. Thus,MP2withasmallbasissetmaybenifit fromthecancelationoferrors. Theauthorsselectedabasis(cc-pVDZ+ 78 )suchthatthe errors from MP2 and the basis cancel out each other. Binding energies of the different benzene dimer configurations at the MP2/cc-pVDZ+ level of theory are within 15% from the best known theoretical values [CCSD(T)/CBS]. This level of theory has been used to investigate the nature of binding in benzene trimers and tetramers. The goal was to understand the importance of the total nonadditive and the individual many- body interaction energies. The authors concluded that the binding energy of trimers and tetramers can be described as a simple nearest neighbor sum with 10% accuracy, and that non-additive (many-body and long range two-body) interactions will become increasinglyimportantforlargerclusters. 38 An important question for the π −π stacked systems is how the share of the non- additive energy component grows with the increase in the number of aromatic rings in the cluster: whether it becomes significant in large clusters or does it stay similar to that of benzene trimers. A related question is if the binding energy of multiple π−π interactions can be well approximated by the sum of the pair-wise nearest neighbor interactionenergies. Toanswerthesefundamentalquestions,oneneedstostudyaseries of the increasingly large aromatic clusters and observe the asymptotic behavior of the non-additive energy component. Unfortunately, aromatic clusters with more than a few benzene rings are beyond the reach of high level ab initio theories. Even using the proposed in Ref. 78 MP2/cc-pVDZ+ method, computing the energy of the cluster with afewtensofthebenzeneringsispresentlyunfeasible. The approach to the multiple π − π interactions problem taken in this work is to employ the effective fragment potential (EFP) — a semi-empirical method designed specifically to account for the inter-molecular interactions. In recent works 27,28,81 , EFP has been used in the extensive studies of benzene dimers. EFP was found to be in excellent agreement with high-level ab initio theories. For the three conformers of the benzene dimer shown in Fig. 1.17, the maximum binding energy was within 0.4 kcal/- mol and the intermolecular separation was within 0.2 ˚ A from the values obtained with CCSD(T)/aug-cc-pVQZ. Additionally, the authors found that the EFP energy profiles (Fig 10.2 and 10.3 in Ref. 28) for the individual types of interactions (electrostatics, polarization, dispersion and exchange-repulsion) along the intermolecular separation coordinateswereinexcellentagreementwiththesymmetry-adaptedperturbationtheory (SAPT).Thus,notonlythetotalbindingenergiesofthebenzenedimersarereproduced wellbyEFP,butalsotheindividualtypesoftheinteractionsareproperlyaccountedfor. 39 ` S T PD R R R 1 R 2 Figure 1.17: Sandwitch (S), t-shaped (T) and parallel-displaced (PD) structures of the benzenedimer. The discussion above indicates that EFP is an excellent candidate for studying the systemswithπ−πinteractions. ThisworkusesournewimplementationofEFPtoinves- tigatemultipleπ−π interactions. Firstly, theEFPbindingenergiesandgeometriesfor thesandwich(S),t-shaped(T)andparallel-displaced(PD)structuresarecomparedwith those obtained at different ab initio levels of theory, including the best known to date theoreticalestimates[CCST(T)/CBS].Then,eightdifferentbenzenetrimersarestudied inordertounderstandhowwellEFPaccountsforthenon-additiveenergycontributions. Finally,theasymptoticbehaviorofmultipleπ−π interactionsisinvestigatedin3types ofbenzeneoligomers. Benzenedimers The parameters for benzene used throughout this work are those recommended in Ref. 28. The geometry of the monomer [d(C-C)=1.3942, d(H-H)=1.0823] computed at the MP2/aug-cc-pVTZ level of theory is from Ref. 82. The effective fragment potential for the benzene monomer has been generated with GAMESS 34 at the HF/6- 311G++(3df,2p) level of theory. As explained in Ref. 26, a basis set of such high qual- ityis necessaryin orderto generateaccuratedispersion parameters. The inputfileused 40 for computing EFP parameters is provided in Appendix A. The parameters have been convertedtotheQ-ChemformatusingthescriptprovidedinAppendixC. The benzene dimer structures shown in Fig. 1.17 where constructed analogously to Ref. 82. The geometries and energies (table 1.1) of 3 benzene dimers corresponding to the maximum binding at the EFP level of theory were found by performing scans alongthedegreesofthefreedomshowninFig.1.17witharesolutionof0.1 ˚ A.Atthese geometries, we verified that our implementation computes the total binding energies andindividualenergycomponentscorrectlybycomparingwithananalogousGAMESS computation. Thegeometriesandmaximumbindingenergieswefoundareonlyslightly different (0.1 ˚ A and 0.15 kcal/mol, respectively) from those obtained in Ref. 28. The discrepancyisduetothedifferenttechniquesoftheelectrostaticsscreening. c Table1.1: Thegeometricparametersandenergiesofthebenzenedimers S T PD Geometricparameters, ˚ A R=4.0 R=5.2 R 1 =3.9,R 2 =1.2 Energy,kcal/mol -1.96 -2.38 -2.18 Table1.2comparestheenergiesofthe3benzenedimerstructuresatdifferentlevels of theory. The best known theoretical estimates are the CCSD(T)/CBS values, which willserveasabenchmarkforourcomparison. WecanseethatCCSD(T)/aug-cc-pVQZ gives the results within 0.15 kcal/mol from the CBS result. Thus, aug-cc-pVQZ can be considered sufficient. However, the MP2 results with this basis set are significantly off (error up to 70%), which is attributed to the overestimation of electron correlation by MP2 in the benzene dimer 83 . The method (MP2/aug-cc-pVDZ+) proposed in Ref. 78 c InRef.28inadditiontothescreeningofcharges,whichourimplementationalsoperforms,thehigher order multipoles (dipoles, quadrupoles, octupoles) were also screened. However, as it was demonstrated inRef.27,theaccuracyofcharge-chargeonlyandcharge-chargeplushigh-orderelectrostaticscreenings issimilar. Therefore,weomittedscreeningofthehigherordermultipolesinourimplementation. 41 is within 0.4 kcal/mol due to the cancellation of errors. Our EFP results are within 0.6 kcal/molandthusareofaccuracycomparabletoMP2/aug-cc-pVDZ+. Table1.2: Dissociationenergiesofthebenzenedimersatdifferentlevelsoftheory. Method a S T PD MP2/aug-cc-pVDZ* -2.83 -3.00 -4.12 MP2/aug-cc-pVTZ -3.25 -3.44 -4.65 MP2/aug-cc-pVQZ* -3.35 -3.48 -4.73 MP2/aug-cc-pVQZ -3.37 -3.54 -4.79 MP2/cc-pVDZ+ -1.87 -2.35 -2.84 CCSD(T)/cc-pVDZ* -1.33 -2.24 -2.22 Est. CCSD(T)/aug-cc-pVQZ* -1.70 -2.61 -2.63 Est. CCSD(T)/CBS -1.81 -2.74 -2.78 EFP -1.96 -2.38 -2.18 a TheresultsfortheabinitiomethodsarequotedfromRef.76 Benzenetrimers Using the geometric parameters obtained for dimers (table 1.1), linear benzene trimers shown in Fig. 1.18 have been constructed. The cyclic structure was obtained by mini- mizingthebindingenergyalongtheintermolecularseparationandrotationintheplane connecting the molecular centers. The minimum corresponds to the molecular sepa- ratioin of 5.3 ˚ A, tilted 15 ◦ from the perpendicular. The geometries of the trimers are providedinAppendixB. Thesetrimerstructuresareanalogoustothetrimersstudiedin Ref.78byMP2. Bindingenergiesforthetrimersalongwiththeirbreakdownintothecomponentsare shownintables1.3and1.4. E 2 istotaltwo-bodyinteractionenergy: E 2 =E 2 (12)+E 2 (13)+E 2 (23), (1.39) 42 Table1.3: Totalandmany-bodyinteractionenergiesofthebenzenetrimers(firstgroup). S PD T1 T2 EFP E 2 (12) -1.96 -2.18 -2.38 -2.40 E 2 (13) 0.01 0.01 0.02 0.00 E 2 (23) -1.96 -2.18 -2.40 -2.38 E 2 -3.91 -4.35 -4.77 -4.79 E 3 0.02 -0.02 0.13 0.06 E total -3.89 -4.38 -4.64 -4.73 E total(dimer) -3.92 -4.36 -4.76 -4.76 E non−add 0.03 -0.02 0.12 0.03 Δ 0.77% -0.37% 2.66% 0.57% MP2/cc-pVDZ+ E 3 0.03 0.00 0.08 0.06 E total -3.83 -5.88 -4.64 -4.73 E total(dimer) -3.74 -5.68 -4.70 -4.70 E non−add -0.09 -0.20 0.06 -0.03 E CP non−add 0.03 -0.06 0.1 0.03 Δ 0.78% -1.02% 2.16% 0.63% CCSD(T)/cc-pVDZ+ E 3 0.04 0.01 0.07 E total -0.90 -1.84 -3.14 E total(dimer) -0.86 -1.72 -3.20 E non−add -0.04 -0.12 0.06 E CP non−add 0.06 0.00 0.10 Δ 6.67% 0.00% 3.18% E 2 (NM) - two body interaction between benzenes N and M, E 2 - total two body interaction, E 3 - three-body interaction, E total - total energy of trimer, E total(dimer) - totalenergyoftrimerestimatedformdimer,E non−add -non-additiveenergycomponent, E CP non−add - non-additive energy component with counterpoise correction, Δ - share of thenon-additiveenergy. 43 Table 1.4: Total and many-body interaction energies of the benzene trimers (second group). C PD-S S-T PD-T EFP E 2 (12) -2.19 -1.96 -1.96 -2.38 E 2 (13) -2.19 0.00 -0.12 -0.11 E 2 (23) -2.19 -2.18 -2.38 -2.18 E 2 -6.57 -4.14 -4.46 -4.67 E 3 -0.24 0.01 -0.07 -0.06 E total -6.81 -4.12 -4.52 -4.73 E total(dimer) -6.57 -4.14 -4.34 -4.56 E non−add -0.25 0.02 -0.18 -0.17 Δ -3.63% 0.43% -4.05% -3.66% MP2/cc-pVDZ+ E 3 -0.33 0.02 -0.03 0.00 E total -7.88 -4.86 -4.49 -5.42 E total(dimer) -7.32 -4.71 -4.22 -5.19 E non−add -0.56 -0.15 -0.27 -0.23 E CP non−add -0.32 -0.01 -0.15 -0.13 Δ -4.06% -0.21% -3.34% -2.40% CCSD(T)/cc-pVDZ+ E 3 -0.25 E total -5.09 E total(dimer) -4.62 E non−add -0.47 E CP non−add -0.26 Δ -5.11% E 2 (NM) - two body interaction between benzenes N and M, E 2 - total two body interaction, E 3 - three-body interaction, E total - total energy of trimer, E total(dimer) - totalenergyoftrimerestimatedformdimer,E non−add -non-additiveenergycomponent, E CP non−add - non-additive energy component with counterpoise correction, Δ - share of thenon-additiveenergy. 44 S T1 T2 PD C PD-T PD-S S-T Figure1.18: Thestructuresofeightbenzenetrimers. whereE 2 (NM)isbindingenergyforthedimers,createdbyremovingonebenzenefrom the trimer and computed in the full basis set of the trimer. E 3 is three-body interaction energycomputedas: E 3 =E total −E 2 (1.40) E total(dimer) is total binding energy estimated as the sum of nearest neighbor pair-wise bindingenergies. E non−add andE CP non−add arethenon-additiveenergycomponentswith- outandwiththecounterpoisecorrection: E non−add =E total −E total(dimer) E CP non−add =E total −E CP total(dimer) (1.41) 45 Δistheshareofthenon-additiveinteractionenergywhichEFPcomputesas: Δ= E non−add E total (1.42) andabinitiomethodsas: Δ= E CP non−add E total (1.43) To evaluate how well EFP describes multiple π−π interactions, it is necessary to conduct a comparison with ab initio methods. The quality of the description of the two-body interactions (between two benzene rings) has been investigated in the exten- sive studies of the benzene dimers energy profiles in Ref. 28. The conclusion was that EFP results are in good agreement with the high level ab initio theories [CCSD(T) and SAPT].Three-bodyinteractionsforEFP,MP2/cc-pVDZ+andCCSD(T)/cc-pVDZ+are compared in tables 1.3 and 1.4. The first important observation is that the ab initio E 3 values quoted from Ref. 78 were computed using the counterpoise correction to account for the basis set superposition error, which results from the basis set incom- pleteness. Three-body interaction energies were calculated by subtracting from total bindingenergyofthetrimerenergiesofthe3possiblebenzenedimerscomputedinthe fullsetbasisofthetrimer. AllEFPmany-bodyinteractionsareduetothepolarizationcomponent. Thus,its E 3 term is three-body polarization energy. From tables 1.3 and 1.4, we can see that three- body interaction energies computed with EFP have the same order of magnitude as for the ab initio methods. For most trimers studied here, the E 3 correction is rather small andmiscellaneouserrorsmayhaveasignificantimpactonitsvalue. Theexceptionisthe cyclic benzene trimer (see Fig. 1.18), which is a true three-body system demonstrating substantial three-body interactions at all the considered levels of theory. For this trimer E 3 computed with EFP is within 0.01 kcal/mol from the CCSD(T) value. The fact that 46 these values almost coincide might be accidental, but their closeness implies a certain levelofaccuracyinthetreatmentofthree-body π−π interactionsbyEFP.Onthebasis of the above comparison, we conclude that EFP provides a good description of three- bodyinteractionsinbenzenetrimers. Anotherinterestingconclusioncanbedrawnfromthecomparisonoftotalnonaddi- tive interaction energies (E non−add ). Firstly, E non−add for EFP are greater than for CCSD(T) and MP2 for all eight considered trimers. However, we concluded that this effect is mostly due to the basis set superposition error (BSSE). E non−add is computed as a difference between trimer binding energy and energies of two dimers (which were derived from the trimer by removing one of the edge benzenes). Thus, the basis sets used for the trimer and the dimers energy computations are different. Since trimer has a larger basis set, its nearest neighborπ−π interactions are stronger than in the corre- sponding dimers. That is, the additional basis functions contribute to the wavefunction relaxation, leading to binding energy lowering. In fact, table 1 of Ref. 78 shows the energies of the dimers computed in the full basis set of the trimer. These energies are 0.05-0.10 kcal/mol lower than for the dimers in its own basis set. Thus, as much as 0.3 kcal/mol for the cyclic trimer and 0.2 kcal/mol for other trimers can be contributed to theloweringofE non−add duetothebasissetrelaxationeffect. SinceEFPdoesnotsuffer fromthebasissetrelaxation,itsE non−add forbenzenetrimersaresystematicallyhigher. E non−add computedbyabinitiomethodscanbecorrectedbyusingthebindingener- gies of the dimers computed in the full basis set of the trimer, similarly to the counter- poisecorrectionoftenemployedtorectifyBSSE.Sincethesevaluesarereadilyavailable intable1ofRef.78,wewereabletocomputecorrectednon-additiveinteractionenergy E CP non−add , which is also provided in tables 1.3 and 1.4. EFP non-additive interaction energies are within 0.03 kcal/mol from E CP non−add for CCSD(T) and within 0.04 kcal/- molfrom MP2values (exceptforthe cyclictrimer). The shareof thetotal non-additive 47 energy correction Δ does not exceed 4% for EFP and MP2. This share is larger for CCSD(T) because the method significantly underestimates the total binding energy. It alsomayattributedtotheinclusionofhigher-orderthree-bodyinteractionenergyterms (e.g.,three-bodydispersion)inCCSD(T). Based on the above comparisons of three-body and total non-additive interaction energies, we conclude that EFP accurately describes the multipleπ−π interactions in theconsideredbenzenetrimers. Largebenzeneoligomers Thelastpartofourinvestigationisastudyoftheasymptoticbehaviorofthenon-additive energy contribution. Using the intermolecular separation parameters for dimers (table 1.1), we constructed the linear oligomers with 10, 20, 30 and 40 benzene molecules. TheoligomersS-N,T-NandPD-N,whereNisthenumberofmonomers,havethesame structural motives as trimers S, T1 and PD shown in Fig. 1.18. The geometries of the oligomers are provided in Appendix B. The total binding energies (E total ), the nearest neighbor estimates (E total(dimer) ), the non-additive energy ( E non−add ) and its share (Δ) inthetotalenergyareshownintable1.5. Thenearestneighborestimatefortheoligomer withNbenzeneringshasbeencomputedas: E total(dimer) =(N −1)E dimer (1.44) From table 1.5, we can see that for all the oligomers E non−add is positive and thus destabilizing in nature, as it was correctly predicted in Ref. 78. It is also in good agreement with the observation that nonadditive energy is positive in the bulk benzene made in Ref. 84. Non-additive energy grows with the size of the system because of the 48 Table1.5: BindingenergiesofbenzeneoligomerscomputedwithEFP. E total E total(dimer) E non−add Δ S-10 -17.14 -17.64 0.50 2.94% S-20 -36.01 -37.24 1.23 3.41% S-30 -54.89 -56.84 1.95 3.56% S-40 -73.76 -76.44 2.68 3.63% T-10 -20.89 -21.42 0.53 2.53% T-20 -44.05 -45.22 1.17 2.65% T-30 -67.21 -69.02 1.81 2.69% T-40 -90.37 -92.82 2.45 2.71% PD-10 -19.54 -19.62 0.08 0.41% PD-20 -41.16 -41.42 0.26 0.63% PD-30 -62.78 -63.22 0.44 0.70% PD-40 -84.40 -85.02 0.62 0.73% a S-N, T-N and PD-N are linear benzene oligomers with N benzene monomers, con- structed using sandwich, t-shaped and parallel-displaced structural motives, respec- tively. largenumberofmany-bodyandlong-rangetwo-bodyinteractionspossibleinthelarger oligomers. The dependence of the share of non-binding energy Δ on the system size is shown in Fig. 1.19. For all 3 types of the oligomers, Δ asymptotically approaches a constant valueanddoesnotexceed4%ofthetotalbindingenergy. Basedontheinvestigationofbenzenetrimersandhigher-orderoligomersattheEFP level of theory, we conclude that the share of total non-additive energy is about 4% forthesystems considered. Thisresultimplies thattotalinteractionenergyin clustered π−πsystemscanbeestimatedreasonablywellasthesumoftwo-bodynearestneighbor interactions. However, it is important to note that neither EFP applied in our work nor MP2 used in the studies of trimers in Ref. 78 account for many-body dispersion. As demonstrated in the recent study 79 , many-body dispersion terms contribute about 49 1 0 2 0 3 0 4 0 0 .5 1 .0 1 .5 2 .0 2 .5 3 .0 3 .5 S h a re o f th e n o n -a d d itiv e e n e rg y , % T h e n u m e r o f b e n z e n e m o n o m e rs S P D T Figure 1.19: Dependence of the non-additive energy share on the length of benzene oligomers. 50% to total three-body interaction energy of the cyclic benzene trimer. Additionally, we expect that the share of the non-additive energy contributions will increase in the three-dimensional clusters and bulk benzene due to the large number of many-body interactions, which can occur between molecules located within a short distance from eachother. 1.8 Conclusions This chapter presented a review of the theoretical approaches to modeling the envi- ronmental effects and a detailed description of the EFP method. The details of our EFP implementation in the Q-Chem package were presented. The developed code was 50 appliedtoinvestigatemultipleπ−π interactionsinbenzeneoligomersfocusingonthe roleofmany-bodyandtotalnon-additiveenergiesinthesesystems. We analyzed previously reported 78 ab initio results for benzene trimers and con- cluded that they underestimate non-additive energy due to the basis set superposition error. Based on the data presented in Ref. 78, we were able to recompute ab initio non-additive energy with the counterpoise correction. Using our EFP code, we com- puted three-body and total non-additive energies for benzene trimers. These energies foralleighttimersincludedinourstudyappearedtobeinexcellentagreementwiththe counterpoise-corrected ab initio values (within 0.04 kcal/mol). The magnitude of non- additiveinteractionenergywaswithin4%oftotalbindingenergyforalltheconsidered benzene clusters. Based on the detailed comparison of the binding energy components (presented in tables 1.3 and 1.4), we concluded that EFP describes the intermolecular π−π interactionswiththeaccuracycomparabletothatofabinitiomethods. Inthefinalpart,ourEFPcodewasappliedtostudylinearbenzeneoligomers. Based ontheresults(presentedintable1.5),weconcludedthattheshareofnon-additiveinter- actionenergyinlinearbenzeneclustersasymptoticallyapproachesaconstantvaluewith the increase in cluster size, and does not exceed 4% of total binding energy. Thus, it is approximatelyofthesamemagnitudeasforthebenzenetrimers. Thisresultimpliesthat therearenosignificantnon-localeffectsinthebindingmechanismofthelinearbenzene clusters. 51 Chapter1References [1] Raghavachari, K. ; Trucks, G.W. ; Pople, J.A. ; Head-Gordon, M. Chem. Phys. Lett.1989,157,479. [2] Bowler, D.R. ; Miyazaki, T. ; Gillan, M.J. J. Phys.: Condens. Matter 2002, 14, 2781. [3] Schwegler,E.;Challacombe,M.J.Chem.Phys.1996,105,2726. [4] White, C.A. ; Johnson, B.G. ; Gill, P.M.W. ; Head-Gordon, M. Chem. Phys. Lett. 1996,253,268. [5] Sch¨ utz,M.;Hetzer,G.;Werner,H.J.Chem.Phys.1999,111,5691. [6] Scuseria,G.E.;Ayala,P.Y.J.Chem.Phys.1999,111,8330. [7] Gogonea,V.;Westerhoff,L.M.;K.M.Merz,Jr.J.Chem.Phys.2000,113,5604. [8] Kitaura, K. ; Ikeo, E. ; Asada, T. ; Nakano, T. ; Uebayasi, M. Chem. Phys. Lett. 1999,313,701. [9] Lippert,G.;Hutter,J.;Parrinello,M.Mol.Phys.1997,92,477. [10] Allunger,N.L.;Yuh,Y.H.;Lii,J.H.J.Am.Chem.Soc.1989,111,8551. [11] Casewit,C.J.;Colwell,K.S.;Rappe,A.K.J.Am.Chem.Soc.1992,114,10035. [12] Mayo,S.L.;Olafson,B.D.;Goddard,W.A.J.Phys.Chem.1990,94,8897. [13] Rizzo,R.C.;Jorgensen,W.L.J.Am.Chem.Soc.1999,121,4827. [14] Firman,T.K.;Landis,C.R.J.Am.Chem.Soc.2001,123,11728. [15] Jorgensen,W.L.;McDonald,N.A.J.Molec.Struct.(Theochem)1998,424,145. [16] MacKerell,A.D.; Bashford,D.; Bellott; Dunbrack,R.L.; Evanseck,J.D.; Field, M.J. ; Fischer, S. ; Gao, J. ; Guo, H. ; Ha, S. ; Joseph-McCarthy, D. ; Kuchnir, L. ; Kuczera, K. ; Lau, F.T.K. ; Mattos, C. ; Michnick, S. ; Ngo, T. ; Nguyen, D.T. ; Prodhom, B. ; Reiher, W.E. ; Roux, B. ; Schlenkrich, M. ; Smith, J.C. ; Stote, R. ; Straub, J. ; Watanabe, M. ; Wiorkiewicz-Kuczera, J. ; Yin, D. ; Karplus, M. J. Phys.Chem.B1998,102,3586. [17] Hill,J.R.;Sauer,J.J.Phys.Chem.1995,99,9536. [18] Lin,H.;Truhlar,D.G.Theor.Chem.Acc.2007,117,185. 52 [19] Gao,J.L.;Truhlar,D.G.Annu.Rev.Phys.Chem.2002,53,467. [20] Sierka,M.;Sauer,J.J.Chem.Phys.2000,112,6983. [21] Monard,G.;Merz,K.M.Acc.Chem.Res.1999,32,904. [22] Day, P.N. ; Jensen, J.H. ; Gordon, M.S. ; Webb, S.P. ; Stevens, W.J. ; M.Krauss; D.Garmer;H.Basch;D.CohenJ.Chem.Phys.1996,105,1968. [23] Jensen,J.H.;Gordon,M.S.J.Chem.Phys.1998,108,4772. [24] Freitag, M.A. ; Gordon, M.S. ; Jensen, J.H. ; Stevens, W.J. J. Chem. Phys. 2000, 112,7300. [25] Gordon, M.S. ; Freitag, M.A. ; Bandyopadhyay, P. ; Jensen, J.H. ; Kairys, V. ; Stevens,W.J.J.Phys.Chem.A2001,105,293. [26] Adamovix,I.;Gordon,M.S.Mol.Phys.2005,103,379. [27] Slipchenko,L.V.;Gordon,M.S.J.Comput.Chem.2007,28,276. [28] Gordon, M.S. ; Slipchenko, L. ; Li, H. ; Jensen, J.H. In Annual Reports in Com- putational Chemistry; Spellmeyer, D.C. , Wheeler, R. , Eds., Vol. 3 of Annual ReportsinComputationalChemistry; Elsevier,2007; pages177–193. [29] Slipchenko,L.V.;Gordon,M.S.Mol.Phys.2009. [30] Li,H.;Netzloff,H.M.;Gordon,M.S.J.Chem.Phys.2006,125. [31] Li,H.;Gordon,M.S.;Jensen,J.H.J.Chem.Phys.2006,124. [32] Li,H.;Gordon,M.S.J.Chem.Phys.2007,126. [33] Li,H.;Gordon,M.S.Theor.Chem.Acc.2006,115,385. [34] Schmidt, M.W. ; Baldridge, K.K. ; J.A. Boatz, S.T. Elbert ; Gordon, M.S. ; J.H. Jensen, S. Koseki ; Mastunaga, N. ; Nguyen, K.A. ; S. Su, T.L. Windus ; Dupuis,M.;Montgomery,J.A.J.Comput.Chem.1993,14,1347. [35] Adamovic,I.;Gordon,M.S.J.Phys.Chem.A2005,109,1629. [36] Balawender,R.;Safi,B.;Geerlings,P.J.Phys.Chem.A2001,105,6703. [37] Chandrakumar, K.R.S. ; Ghanty, T.K. ; Ghosh, S.K. ; Mukherjee, T. J. Molec. Struct.(Theochem)2007,807,93. [38] Day,P.N.;Pachter,R.J.Chem.Phys.1997,107,2990. 53 [39] Merrill,G.N.;Webb,S.P.J.Phys.Chem.A2004,108,833. [40] Yoo,S.;Zahariev,F.;Sok,S.;Gordon,M.S.J.Chem.Phys.2008,129,8. [41] Freitag,M.A.;Gordon,M.S.AbstractsofPapersoftheAmericanChemicalSoci- ety1999,218,U511; 1. [42] Merrill,G.N.;Gordon,M.S.J.Phys.Chem.A1998,102,2650. [43] Merrill,G.N.;Webb,S.P.J.Phys.Chem.A2003,107,7852. [44] Merrill,G.N.;Webb,S.P.;Bivin,D.B.J.Phys.Chem.A2003,107,386. [45] Netzloff,H.M.;Gordon,M.S.J.Chem.Phys.2004,121,2711. [46] Day,P.N.;Pachter,R.andGordon,M.S.;Merrill,G.N.J.Chem.Phys.2000,112, 2063. [47] Shao, Y. ; Molnar, L.F. ; Jung, Y. ; Kussmann, J. ; Ochsenfeld, C. ; Brown, S. ; Gilbert, A.T.B. ; Slipchenko, L.V. ; Levchenko, S.V. ; O’Neil, D. P. ; Distasio Jr., R.A.;Lochan,R.C.;Wang,T.;Beran,G.J.O.;Besley,N.A.;Herbert,J.M.;Lin, C.Y. ; Van Voorhis, T. ; Chien, S.H. ; Sodt, A. ; Steele, R.P. ; Rassolov, V. A. ; Maslen,P.;Korambath,P.P.;Adamson,R.D.;Austin,B.;Baker,J.;Bird,E.F.C. ; Daschel, H. ; Doerksen, R.J. ; Drew, A. ; Dunietz, B.D. ; Dutoi, A.D. ; Furlani, T.R. ; Gwaltney, S.R. ; Heyden, A. ; Hirata, S. ; Hsu, C.-P. ; Kedziora, G.S. ; Khalliulin, R.Z. ; Klunziger, P. ; Lee, A.M. ; Liang, W.Z. ; Lotan, I. ; Nair, N. ; Peters, B. ; Proynov, E.I. ; Pieniazek, P.A. ; Rhee, Y.M. ; Ritchie, J. ; Rosta, E. ; Sherrill,C.D.;Simmonett,A.C.;Subotnik,J.E.;WoodcockIII,H.L.;Zhang,W. ;Bell,A.T.;Chakraborty,A.K.;Chipman,D.M.;Keil,F.J.;Warshel,A.;Herhe, W.J.;SchaeferIII,H.F.;Kong,J.;Krylov,A.I.;Gill,P.M.W.;Head-Gordon,M. Phys.Chem.Chem.Phys.2006,8,3172. [48] Saenger, W. Principles of Nucleic Acid Structure; Springer-Verlag: New York, 1984. [49] Buckingham,A.D.Q.Rev.Chem.Soc.1959,13,183. [50] Stone,A.J.Chem.Phys.Lett.1981,83,233. [51] Mulliken,R.S.J.Chem.Phys.1955,23,1833. [52] Stone,A.J.;Alderton,M.Mol.Phys.1985,56,1047. [53] Stone, A.J. In The Theory of Intermolecular Forces; Oxford University Press, 1997; pages50–62. [54] Stone,A.J.J.Chem.TheoryComput.2005,1,1128. 54 [55] Rabinowitz, J.R. ; Namboodiri, K. ; Weinstein, H. Int. J. Quant. Chem. 1986, 29, 1697. [56] Wolfsberg,M.;Helmholz,L.J.Chem.Phys.1952,20,837. [57] Boys,S.F.Rev.Mod.Phys.1960,32,296. [58] Boys,S.F.InQuantumScienceofAtoms,Molecules,andSolids; AcademicPress, 1966; pages253–262. [59] Sanz-Garcia,A.;Wheatley,R.J. ChemPhysChem2003,5,801. [60] Tang,K.T.;Toennies,J.P.J.Chem.Phys.1984,80,3726. [61] Jensen,J.H.;Gordon,M.S.Mol.Phys.1996,89,1313. [62] ISO/IEC ISO/IEC 14882:2003: Programming languages: C++ Technical report, InternationalOrganizationforStandardization,Geneva,Switzerland.,2003. [63] Press,W.H.;Teukolsky,S.A.;Vetterling,W.T.;Flannery,B.P.NumericalRecipes 3rd Edition: The Art of Scientific Computing; Cambridge University Press: New York,NY,USA,2007. [64] Rossum, G.V. The Python Language Reference Manual; Network Theory Ltd., 2003. [65] Hunter,C.A.;Sanders,J.K.M.J.Am.Chem.Soc.1990,112,5525. [66] Burley,S.K.;Petsko,G.A.Science1985,229,23. [67] Brana, M.F. ; Cacho, M. ; Gradillas, A. ; Pascual-Teresa, B. ; Ramos, A. Curr. Pharm.Des.2001,7,1745. [68] Gazit,E.Chem.Soc.Rev.2007,36,1263. [69] Desiraju,G.R.;Gavezzotti,A.JCSChem.Comm.1989,page621. [70] Askew,B.;Ballester,P.;Buhr,C.;Jeong,K.S.;Jones,S.;Parris,K.;Williams, K.;Rebek,J.J.Am.Chem.Soc.1989,111,1082. [71] Smithrud,D.B.;Diederich,F.J.Am.Chem.Soc.1990,112,339. [72] Muller-Dethlefs,K.;Hobza,P. Chem.Rev.2000,100,143. [73] Felker,P.M.;Maxton,P.M.;Schaeffer,M.W.Chem.Rev.1994,94,1787. [74] Leopold, K.R. ; Fraser, G.T. ; Novick, S.E. ; Klemperer, W. Chem. Rev.1994, 94, 1807. 55 [75] Mueller-Dethlefs,K.;Dopfer,O.;Wright,T.G. Chem.Rev.1994,94,1845. [76] Sinnokrot,M.O.;Sherrill,C.D.J.Phys.Chem.A2006,110,10656. [77] Ye,X.;Li,Z.;Wang,W.;Fan,K.;Xu,W.;Hua,Z.Chem.Phys.Lett.2004,397, 56. [78] Tauer,T.P.;Sherrill,C.D.J.Phys.Chem.A2005,109,10475. [79] Podeszwa,R.;Szalewicz,K.J.Chem.Phys.2007,126,194101. [80] Podeszwa,R.J.Phys.Chem.A2008,112,8884. [81] Smith,T.;Slipchenko,L.V.;Gordon,M.S.J.Phys.Chem.A2008,112,5286. [82] Sinnokrot,M.O.;Valeev,E.F.;Sherrill,C.D.J.Am.Chem.Soc.2002,124,10887. [83] Tsuzuki, S. ; Uchimaru, T. ; Matsumura, K. ; Mikami, M. ; Tanabe, K. Chem. Phys.Lett.2000,319,547. [84] Podeszwa,R.;Bukowski,R.;Szalewicz,K.J.Phys.Chem.A2006,110,10345. 56 Chapter2 Structure,vibrationalfrequencies, ionizationenergies,andphotoelectron spectrumofthepara-benzyneradical anion 2.1 Introduction Para-benzyne 1 , an intermediate in Bergman cyclization 2 , is believed to be a warhead in antitumorenediyneantibiotics 3–5 duetoitsabilitytocausedouble-strandDNAcleavage andself-programmedcelldeathorapoptosis 6,7 ,andextensivestudieshavebeencarried outinthepast20yearstounderstandthefactorsthatcontrolthereaction 8–10 . Originally, it was interest in para-benzyne electronic structure that motivated the development of procedures to generate the corresponding anion for use in photoelectron spectroscopy experiments 11,12 . Recent studies suggested that anionic cycloaromatization reactions arealsofeasible 13,14 includingtheanionicversionoftheBergmancyclization 15 . Larger sensitivity to substituents 15 suggests that the reactions of the anions might be designed tobemoreselectivethantheprototypicnon-anionicBergmancyclization,importantfor theirapplicationsincancer-DNAdamagingdrugs. Previous theoretical studies of the para-benzyne anion 16 predicted strong vibronic interactions(oftenreferredtoaspseudoorsecond-orderJahn-Tellereffect)thatproduce 57 alowersymmetryC 2v structureinadditiontotheD 2h minimum. Theenergydifference between the two structures was found to be highly sensitive to the correlation method, with higher-level models favoring D 2h . The analysis of the photoelectron spectrum of the anion 12 supported the assignment of the D 2h symmetry, and it was concluded that the lower symmetry structures are artifacts of incomplete correlation treatment, a well known phenomenon in electronic structure 17 . Interestingly, density functional calcula- tions produced only D 2h structures, which was regarded as success of DFT methodol- ogy and the authors concluded that “this level of theory is particularly well-suited for computational studies of distonic radical anions derived from diradicals” 16 . The robust behavior of DFT in cases of symmetry breaking has been noted by other researches as well 18 . The goal of this work is to investigate the performance of coupled-cluster (CC) and equation-of-motion(EOM)methodsinthischallengingcaseofstrongvibronicinterac- tions and to assess the quality of the ab initio and DFT results by comparing the cal- culated photoelectron spectrum with the experimental one. We quantify the symmetry loweringbytheselectedgeometricalparameters,chargedistributions,andrelativeener- gies. Inaddition,theshapeofthepotentialenergysurfacealongthesymmetrybreaking coordinate is analyzed. We also present an efficient scheme for calculating ionization energies(IEs)oftheanioninthespiritofisodesmicreactions. Severalrecentstudies 19–21 discusseddifferentaspectsofsymmetrybreakinginelec- tronic structure calculations distinguishing between: (i) purely artefactual spatial and spin symmetry breaking of approximate wave functions, i.e., L¨ owdin dilemma 22 ; (ii) real interactions between closely-lying electronic states that result in lower-symmetry structures, significant changes in vibrational frequencies (see Fig. 1 from Ref. 21), and evensingularities(firstorderpoles)inforceconstants(e.g.,Fig. 3fromRef.20). While thelatterisarealphysicalphenomenon,itmayormaynotbeaccuratelydescribedbyan 58 approximate method. If vibronic interactions are overestimated, calculations may yield incorrectlower-symmetrystructure,andviceverse. Formally, the effect of the interactions between electronic states on the shape of potential energy surface can be described 20 by the Herzberg-Teller expansion of the potentialenergyofelectronicstatei: V i =V 0 + X α <Ψ i | ∂H ∂Q α |Ψ i >Q α + 1 2 X α <Ψ i | ∂ 2 H ∂ 2 Q α |Ψ i >Q 2 α − X α X j6=i <Ψ i | ∂H ∂Qα |Ψ j > 2 E j −E i Q 2 α , (2.1) where Q α denotes normal vibrational modes, wave functions Ψ k and energies E k are adiabatic wave functions and electronic energies at Q α = 0. The last term, which is quadraticinnucleardisplacement—hencesecond-orderJahn-Teller,describesvibronic effects. Forthegroundstate,i.e.,whenE j >E i ,itcausessofteningoftheforceconstant along Q α , which may ultimately result in a lower symmetry structure, e.g., see Fig. 1 from Ref. 21. The above mentioned poles in force constants originate in the energy denominator. The formal analysis of the CC and EOM-CC second derivatives by Stanton 20 demonstrated that the quadratic force constant of standard coupled-cluster methods, e.g., coupled-cluster with singles and doubles (CCSD), necessarily contains unphysi- cal terms, which may become significant and spoil the CC force constant in the cases of strongly interacting states. Nevertheless, the standard CC methods are quite reliable for systems exhibiting pseudo Jahn-Teller interactions, as follows from Refs. 20,21, as well as a host of numerical evidence. EOM-CC methods, on the other hand, provide mostsatisfactorydescriptionoftheinteractingstates 20 ,notatallsurprisinginviewofa 59 multistage nature of EOM. It should be noted that both CC and EOM-CC may exhibit spurious frequencies if the orbital response poles, i.e., coupled-perturbed Hartree-Fock (CPHF) poles, are present due to near-instabilities in the reference, a condition closely related to (i). Perturbative corrections, e.g., as in CCSD(T), do not provide systematic improvement,andoftenresultsinwiderinstabilityvolcanoes 19 . By analyzing the properties of linear response of multi-configurational SCF (MCSCF)andtherelationshipbetweenresponseequationsandpoles’structure,Stanton dispelled a common misconception of relying on multi-reference treatment of systems with strong vibronic interactions 20 . Using the same arguments, he also concluded that DFT describes pseudo Jahn-Teller effects reasonably well. While some numerical evi- dence 18 seemstosupportthisconclusion,otherexamples 21 indicatedthatDFT(B3LYP) tendtounderestimatethemagnitudeofvibronicinteractions. Ourresultsareinfavorof thelatter,andweattributethisbehaviortoself-interactionerror. The structure of this chapter is as follows. The next section describes molecular orbitals (MOs) and relevant electronic states of the para-benzyne anion, as well as the natureofvibronicinteractionsinthissystem. Section2.3outlinescomputationaldetails. Section2.4discussesthecompetitionbetweenlowerandhighersymmetrystructuresat differentlevelsoftheory. Photodetachmentspectraarepresentedinsection2.5. Differ- ent computational strategies for calculating IEs are discussed in section 2.7. Our final remarksaregiveninthelastsection. 2.2 Molecularorbitalpictureandtheoriginofvibronic interactionsinthepara-benzyneanion ThefrontierMOsofthepara-benzyneanionandthediradicalareshowninFig.2.1. At D 2h , the two orbitals belong to different irreps, b 1u anda g , and their density is equally 60 distributed between the two radical centers (C 1 and C 4 ). At C 2v , both orbitals become a 1 , and their densities are localized at different carbons. As in the benzyne diradical, thelowestMO,whichhoststwoelectronsintheanion,isofanti-bondingcharacterwith respect to C 1 -C 4 , whereas the singly occupied MO is bonding. In the case of the dirad- ical, three singlet states and one triplet state are derived from different distributions of twoelectronsonthesetwoorbitals: X 1 A g ,2 1 A g ,1 1 B 1u ,2 3 B 1u . Thegroundsingletstate belongstothesameirrepasthedoublyexcitedsinglet,andbothhavesignificantmulti- configurationalcharacter. Likewise,distributingthreeelectronsonthesetwoorbitalsin theaniongivesrisetotwoelectronicconfigurationsshowninthelowerpanelofFig.2.1. AtD 2h ,thesetwodeterminantsareofdifferentsymmetryandcorrespondtotwodistinct electronic states, X 2 A g and 2 B 1u , separated by the energy gap of about 0.95 eV, as cal- culatedbyEOM-CCforelectron-attachedstates(seebelow). At C 2v ,bothdeterminants are of A 1 symmetry and, therefore, can interact. This interaction increases the energy separationbetweenthetwoelectronicstatesthatacquiremulti-configurationalcharacter and,therefore,“softens”thecorrespondingnormalmodeinthelowerstate. Iftheinter- actionissufficientlystrong,thatis,ifthecouplingmatrixelementislargerelativetothe energygapbetweentheinteractingstates,thelowerpotentialenergysurface(PES)may developalowersymmetryminimum. Inotherwords, loweringsymmetrystabilizesthe lowest electronic state through configurational interaction, and this is the driving force forsymmetrybreaking. This analysis of the wave functions of interacting states complements the lessons learnedfromEq.(2.1)andpoles’structureofquarticforceconstants 20,21 . Forexample, it shows that two factors are important for accurate description of systems with consid- erable vibronic interactions: (i) balanced description of possibly multi-configurational wave function (non-dynamical correlation); and (ii) quantitative accuracy in describing 61 C 2v 2h D a 1 a 1 a g b 1u C1 C1 C1 C1 C4 C4 C4 C4 X 2 A g A 2 B 1u X 2 A 1 A 2 A 1 6a g 5b 1u 6 5 6a g 5b 1u g 5 11a 1 10a 1 11 10 11a 1 10a 1 11 10 Figure 2.1: Frontier MOs and leading electronic configurations of the two lowest elec- tronicstatesatD 2h andC 2v structures. 62 excited state ordering (dynamical correlation). Unbalanced treatment may overempha- size the magnitude of vibronic interaction and, therefore, may lead to so-called arte- factual symmetry breaking 17 . While (i) can be achieved by multi-configurational self- consistent field (MCSCF), (ii) is much more difficult to satisfy. For example, MCSCF exhibits artificial symmetry breaking in NO 3 , which is removed when dynamical cor- relation is included 23 . Ground-state coupled-cluster methods, e.g., CCSD or CCSD(T), may fail due to the insufficient description of non-dynamical correlation, as it happens inNO 3 24,25 ,althoughtheyoftentacklesuchchallengingsystemssuccessfully 26 . The analysis of the MOs and the electronic states of the para-benzyne anion (Fig. 2.1) suggests a moderate multi-configurational character of the ground state wave function, and this is confirmed by EOM-EA amplitudes (as calculated at EOM-EA- CCSD/6-311+G**) For the D 2h structure, the EOM coefficient of the leading config- uration (the A g deteminant from Fig. 2.1 and 2.2) is 0.86. The three next important configurations have weights of 0.28, 0.26 and 0.07. Likewise, for theC 2v structure the weights of the four leading configurations are 0.87, 0.27, 0.16 and 0.15. The EOM- SF-CCSD/6-311+G**amplitudesaresimilar: theweightsoftheleadingconfigurations are 0.81 and 0.80 for the D 2h and C 2v structures, respectevily. This confirms minor multi-configurationalcharacterofthepara-benzyneanion. Forcomparison,theweights ofthefourleadingconfigurationsofthesingletpara-benzynediradical,whichisatruly multiconfigurational system, are 0.57, 0.38, 0.32 and 0.31 and the two leaing configu- rations are double excitations with respect to each other. Thus, we expect ground-state CC methods to be capable of treating C 6 H − 4 reasonably well, unlike heavily multicon- figurationalwavefunctions,asthoseofthesingletdiradicals 24 . The EOM methods that simultaneously include both dynamical and non-dynamical correlation are particularly attractive for the systems with vibronic interactions. More- over, as formal analysisof Stanton shows, the pole structure of the exact force constant 63 6a g 5b 1u 6 5 6a g 5b 1u g 5 6a g 5b 1u g 5 6a g 5b 1u g 5 Reference: 7a g g p-C H 6 2- 4 Quartet Triplet p-C H 6 4 p-C H 6 - 4 Singlet + Doublet p-C H 6 - 4 Ψ target = R IP Ψ ref Ψ target = R EA Ψ ref Ψ target = R SF Ψ ref Figure 2.2: The two lowest electronic states of C 6 H − 4 can be described by: (i) EOM-IP with the dianion reference; (ii) EOM-EA with the triplet neutral p-benzyne reference; (iii)EOM-SFwiththequartetp-benzyneanionreference. of Eq. (2.1) is correctly reproduced within EOM-CC formalism. From the configura- tion interaction point of view, the balanced description of two important determinants from Fig. 2.1 is easily provided by several EOM-CC methods, i.e., ionization poten- tial, electron-attachment, and spin-flip EOM-CC (EOM-EA, EOM-IP, and EOM-SF, respectively),asexplainedinFig.2.2. AllthreeEOMmodelsincludedynamicalcorre- lationaswell. Theirreliabilitydependsonhowwellthecorrespondingreferencestateis described by the single-reference CCSD method and possible (near)-instabilities of the Hartree-Fockreferences. 64 2.3 Computationaldetails Two basis sets were employed in this work. All the geometry optimizations, fre- quencies’ calculations and most single point calculations were performed with the 6- 311+G** 27,28 basis. AdditionalsinglepointcalculationsofIEsusedtheaug-cc-pVTZ 29 basis. The choice of reference in our calculations requires additional comments. The ground-statemethods[i.e.,MP2,CCSD,CCSD(T)],whichemploy 2 A g reference,were conducted using pure-symmetry ROHF and UHF references. The < S 2 > values for UHFwere0.869and1.226atD 2h andC 2v ,respectively. Thestabilityanalysisrevealed severalsymmetry-breakingtripletinstabilitiesat D 2h . Thelowest-energyUHFsolution is heavily spin-contaminated, < S 2 >=1.248, and is also of broken symmetry. When employedintheCCSDcalculations,thelower-energybroken-symmetryUHFreference yielded higher correlated energies, in agreement with our experience (see also foot- note 17 in Ref. 20). Thus, we consider symmetry-pure UHF references to be more appropriateinCCcalculations,eventhoughthisresultsincuspsinthePESscansalong the D 2h → C 2v distortion. Using symmetry-broken reference at D 2h yields continuous curves. In the cases of moderate spin-contamination (as for pure-symmetry UHF refer- ences), CCSD is orbital insensitive, however, the (T) correction is often more reliable for ROHF (see, for example, results from Ref. 30). The optimized-orbital CC models avoidambiguityofthereferencechoicefornon-variationalCCSD,and,inviewofUHF instabilities, we consider optimized orbitals coupled-cluster with doubles (OO-CCD) referencetobemostappropriate(amongtheground-stateCCmethods)forthissystem. References for the EOM calculations (see Fig. 2.2) were not well behaved either. Theclosed-shell 1 A 1 referenceforEOM-IPcorrespondstotheC 6 H 2− 4 dianion,whichis unstable. This caused severe convergence problems and we abandoned EOM-IP treat- ment. TheEOM-EAcalculationsemployedthe 3 B 1u state( 3 A 1 atC 2v )ofthediradical. 65 Unexpectedly, the corresponding UHF reference develops strong spin contamination at relatively small displacements from D 2h : for example, < S 2 >=2.026 and 2.468 andD 2h andC 2v , respectively. This compromises the reliability of EOM-EA-CCSD at C 2v . The problem can be resolved by performing the EOM-EA calculations using the OO-CCDtripletreference(EOM-EA-OD).TheEOM-SFcalculationsemploythe 4 B 1u reference( 4 A 1 atC 2v ),whichcanbederivedfromthe 3 B 1u configurationbyplacingan additionalelectrononaσ-like a g orbital. Despitethesimilaritytotheparenttriplet,the quartet reference has been found to be well behaved throughout the whole range of the D 2h →C 2v distortions,e.g.,<S 2 >=3.776and3.779andD 2h andC 2v ,respectively. Optimized geometries were obtained at the UHF, ROHF, UHF-CCSD 31 , UHF- CCSD(T) 32 ,OO-CCD 31,33 ,EOM-EA-CCSD 34 ,EOM-SF-CCSD 35,36 ,andB3LYP 37 lev- els of theory. The latter two methods have produced only D 2h structures. For all other methods, two minima were found, and the C 2v - D 2h energy differences ΔE’s were cal- culated from the corresponding total energies. In addition, ΔE’s were calculated by a variety of methods at the UHF-CCSD optimized geometries. The lowest electron- ically excited states for C 2v and D 2h structures where computed by EOM-EA-CCSD at the ground state geometries optimized at the same method. We also performed theUHF-CCSD,UHF-CCSD(T),OO-CCD,OO-CCD(T),EOM-EA-CCSD,EOM-SF- CCSD and EOM-EA-OD PES scans along the D 2h - C 2v displacement coordinateR(x) definedas: R(x)=(1−x)∗r(D 2h )+x∗r(C 2v ), (2.2) wherex is a fraction ofC 2v structure and vectorr denotes Cartesian coordinates of the D 2h andC 2v structures. The photodetachment spectra were calculated within double harmonic and parallel normalmodesapproximations(usingnormalmodesofanion)byprogramPES 38 .Inall spectracalculations, weemployedtheSF-DFTgeometryandfrequenciesofthe singlet 66 stateofthediradical 39 ,whilethetripletstatewasdescribedbyCCSD.Fortheanion,we considered UHF-CCSD and B3LYP geometries and harmonic vibrational frequencies oftheanion. Charge distributions were obtained from the natural bond orbital (NBO) analysis 40 . Dipolemomentswerecomputedrelativetothemolecularcenterofmass. TheMulliken conventions 41 for symmetry labels were used: for both C 2v and D 2h , the plane of the molecule is YZ and the Z-axis passes through the two carbons that host the unpaired electrons. ElectronicstructurecalculationswereperformedusingQ-Chem 42 andACES II 43 . 2.4 Equilibrium structures, charge distributions, and the D 2h → C 2v potential energy profiles at different levelsoftheory This section discusses the D 2h and C 2v optimized structures, respective energy differ- ences, PES scans, as well as changes in charge distributions upon symmetry lowering. AspointedoutbyNashandSquires 16 ,mosttheoreticalmethodspredicttwominimaon the PES, with D 2h and C 2v symmetries. As follows from the MOs (see Fig. 2.1), the latterstructuresarecharacterizedbyalocalizedcharge,i.e.,areofadistonictype. The optimized D 2h structures are presented in Fig. 2.3. Fig. 2.4 and Fig. 2.5 sum- marizes differences between theD 2h andC 2v structures calculated at different levels of theory. Forcoupled-clustermethods,thechangesingeometriesbetweendifferentlevels of theory (Fig. 2.4, 2.5) are not significant and are smaller for the D 2h structure than for C 2v . For example, the variation in bond lengths and angles for D 2h structures are 67 H H H H 1.391 1.398 1.384 1.390 1.393 1.393 1.093 1.096 1.092 1.093 1.093 1.092 121.9 121.8 122.0 121.9 121.6 121.5 1.435 1.440 1.433 1.435 1.433 1.433 116.5 116.5 117.2 116.8 116.6 116.6 CCSD CCSD(T) B3LYP OO-CCD EA-CCSD SF-CCSD Figure2.3: D 2h structuresoptimizedatdifferentlevelsoftheory. below0.02 ˚ Aand2 ◦ ,respectively. ThedifferencesarelargerfortheC 2v structures,pos- sibly because of larger spin-contamination of the doublet and triplet UHF references. TheD 2h structurescalculatedbytheCCmethodswithdoublesubstitutions,i.e.,CCSD, OO-CCD, EOM-EA-CCSD, and EOM-SF-CCSD, are very close to each other, which isreassuringinviewoftheinstabilityofthedoubletUHFreference. Thechangesinthe bondlengthsandanglesbetweenCCSD(T)andCCSDarewithin0.01 ˚ Aand1 ◦ ,which is consistent with typical differences for well behaved molecules (Ref. 44). The largest change due to (T) is observed for the distance between the radical centers. The B3LYP structure is considerably different, e.g., the corresponding C 1 -C 4 distance is 0.035 ˚ A shorter than the CCSD(T) value. The above trends suggest that the accuracy of the CC structures is not affected by the HF instabilities and vibronic interactions, and question thereliabilityofB3LYPresults. 68 B 3 L Y P C C S D O O -C C D C C S D (T ) E A -C C S D S F -C C S D 1 1 2 1 1 4 1 1 6 1 1 8 1 2 0 1 2 2 A n g le , d e g M e th o d C 2 v : C 4 D 2 h : C 1 C 2 v : C 1 B 3 L Y P C C S D O O -C C D C C S D (T ) E A -C C S D S F -C C S D 1 .3 8 1 .3 9 1 .4 0 1 .4 1 1 .4 2 1 .4 3 1 .4 4 1 .4 5 B o n d le n g th , A M e th o d D 2 h : C 2 -C 3 C 2 v : C 2 -C 3 C 2 v :C 1 -C 2 D 2 h : C 1 -C 2 C 2 v : C 3 -C 4 Figure 2.4: CCC angles, bonds lengthscalculated by different coupled-cluster methods andB3LYP. The comparisons between the D 2h and C 2v structures allows one to quantify the degree of symmetry lowering at each level of theory. For example, the magnitude of 69 B 3 L Y P C C S D O O -C C D S F -C C S D E A -C C S D C C S D (T ) 2 .8 7 2 .8 8 2 .8 9 2 .9 0 2 .9 1 D is ta n c e , A M e th o d D 2 h : C 1 -C 4 C 2 v : C 1 -C 4 Figure 2.5: The distance between the radical-anion centers calculated by different coupled-clustermethodsandB3LYP. C 2v distortionscanbecharacterizedbythedifferencebetweenC 1 -C 2 andC 3 -C 4 ,which are identical in D 2h . As shown in the lower panel of Fig. 2.4, the difference in C 1 -C 2 between theD 2h andC 2v structures is largest for CCSD and OO-CCSD and is reduced bytheinclusionoftriples. AttheEOM-EA-CCSDlevel,itshrinksevenfurther. TheC 1 - C 2 —C 2 -C 3 difference,aswellasdifferenceintheC 1 andC 4 angles,followsthesame trend. Thus, the magnitude of symmetry breaking is smaller for higher-level methods, or methods that use better references. This, along with the absence of a C 2v minimum attheEOM-SF-CCSDlevel,suggeststhatitsoccurrenceisartefactual. The degree of symmetry breaking can also be quantified by the charge localization, which gives rise to a non-zero dipole moment shown in Fig. 2.6. The correspond- ing atomic charges calculated using the NBO procedure are consistent with the dipole moment changes, except for ROHF, which may be an artifact of the NBO analysis. In agreement with the trends in structural changes, the charge localization is smaller for 70 -0 .4 8 -0 .4 4 -0 .4 0 -0 .3 6 -0 .3 2 R O H F -E A -C C S D q (C 1 ), a .u U H F R O H F U H F -C C S D R O H F -C C S D U H F -E A -C C S D 2 .0 2 .5 3 .0 3 .5 4 .0 4 .5 5 .0 5 .5 D ip o le , d e b y e U H F R O H F U H F -C C S D R O H F -C C S D U H F -E A -C C S D R O H F -E A -C C S D Figure 2.6: The NBO charge at C 1 (top panel) and dipole moments (bottom panel) at the optimized C 2v geometries. The ROHF-CCSD and ROHF-EA-CCSD values are calculatedatthegeometriesoptimizedbythecorrespondingUHF-basedmethod. 71 more correlated wave functions. Interestingly, the overall dipole moment decreases in spite of increased C 1 -C 4 distance. Thus, the charge localization patterns also suggest thatthesymmetryloweringisanartifactofanincompletetreatmentofelectroncorrela- tion. Finally, symmetry lowering can be characterized by energy differences (ΔE’s) between theD 2h andC 2v structures summarized in Table 2.1 and Fig. 2.7. In addition, severalPESscansareshowninFig.2.8. Thefirstpartofthetableandtheupperpanelof Fig.2.7summarizesinglepointcalculationsattheCCSDoptimizedstructures(MCSCF andCASPT2NemploytheMCSCFgeometriesandarefromRef.16),whilethesecond part and the lower panel of Fig. 2.7 present energy differences calculated using the respective optimized structures. All ΔE’s from Table 2.1 and Fig. 2.7 are calculated using symmetry-pure UHF references at D 2h (see section 2.3). Fig. 2.8 shows both the symmetry-pure and symmetry-broken CCSD/CCSD(T) values at D 2h . For these meth- ods, the correct-symmetry reference yields lower correlated energies, which results in cusps on the PESs (see section 2.3). OO-CCD has only one solution at D 2h minimum, andtherespectivecurvesaresmooth. Similarly to the trends in structures and charge localization, the energy difference between the D 2h and C 2v structures is also very sensitive to correlation treatment, in agreement with earlier studies 16 . UHF and ROHF place the C 2v structure significantly lower (by about 1 eV), than the D 2h one, while MP2 gives the almost exactly oppo- site result. At higher level methods, the energy difference becomes much smaller. At the CCSD and OO-CCD levels, the C 2v structure is just a little (<0.1 eV) more stable than the D 2h one. Finally, (T) correlation at the CCSD or OO-CCD levels changes the balanceandtheD 2h structurebecomestheglobalminimum. Similarly,attheEOM-EA- CCSD level, D 2h is also the lowest minimum. The B3LYP and EOM-SF PES do not haveaC 2v minimum. Overall,thePESalongthisdistortionisratherflat(duetovibronic 72 -1 .2 -0 .8 -0 .4 0 .0 0 .4 0 .8 1 .2 U H F -E A -O D O O -C C D O O -C C D (T ) R O H F -C C S D (T ) E (C 2 v ) - E (D 2 h ), e V U H F R O H F U H F -M P 2 R O H F -M P 2 U H F -C C S D R O H F -C C S D D F T /B 3 L Y P U H F -C C S D (T ) U H F -E A -C C S D U H F -S F -C C S D M C S C F C A S P T 2 N -1 .2 -0 .8 -0 .4 0 .0 0 .4 0 .8 1 .2 O O -C C D E (C 2 v ) - E (D 2 h ), e V U H F R O H F U H F -C C S D U H F -C C S D (T ) U H F -E A -C C S D Figure 2.7: Energy differences between theC 2v andD 2h structures (6-311+G** basis). Upper panel: calculated using the UHF-CCSD optimized geometries. MCSCF and CASPT2N employ the MCSCF geometries. Lower panel: calculated at the respective optimizedgeometries. 73 Table2.1: EnergydifferencesbetweentheC 2v andD 2h minima(eV)inthe6-311+G** basis set. Negative values correspond to the C 2v minimum being lower. ZPEs are not included. Method UHF ROHF atUHF-CCSDoptimizedgeometries HF -0.999 -1.049 MP2 1.191 0.887 CCSD a -0.046(-0.116) -0.130 CCSD(T) a 0.140(-0.052) 0.077 B3LYP 0.070 5050 b -0.149 OO-CCD -0.096 OO-CCD(T) 0.081 EOM-EA-CCSD 0.086 EOM-EA-OD 0.050 EOM-SF-CCSD 0.052 atthegeometriesoptimizedbythecorrespondingmethod HF -0.971 -1.039 5050 -0.136 CCSD -0.046 CCSD(T) 0.140 OO-CCD -0.097 EOM-EA-CCSD 0.071 a Calculatedusingsymmetry-pureandsymmetry-broken(valuesinparenthesis)doublet UHFreferencesatD 2h . b DFTfunctionalwith50%oftheHFexchange. interactions): themostreliableestimatesofΔEsbetweenthetwoCCSDoptimizedmin- imarangebetween0.05eV(EOM-SF-CCSD)to0.08eV[CCSD(T)andOO-CCD(T)]. Optimized-orbitalandregularCCmethodsgivesimilarenergyorderings,althoughOO- CCDandOO-CCD(T)placethe D 2h structurealittlehigherthanthecorrespondingCC methodsbasedontheHartree-Fockreference. EOM-SFbehavesimilarlytoEOM-EA. 74 0 2 0 4 0 6 0 8 0 1 0 0 -0 .1 6 -0 .1 2 -0 .0 8 -0 .0 4 0 .0 0 0 .0 4 0 .0 8 0 .1 2 E -E (C 2 v ), e V % o f C 2 v C C S D C C S D (T ) O O -C C D O O -C C D (T ) E A -C C S D S F -C C S D E A -O D B 3 L Y P Figure 2.8: The potential energy profiles along the D 2h - C 2v scans. At the D 2h struc- ture, CCSD and CCSD(T) energies are calculated using symmetry-broken and pure- symmetryUHFreferences(seesection2.3). Similar trends, although with larger variations, were observed for multi-reference methods, e.g. MCSCF and CASPT2. The previous study 16 reported the energy dif- ferences between the C 2v and D 2h structures calculated by multi-reference methods as -0.286 eV for MCSCF(9,8)/cc-pVDZ and 0.388 for CASPT2N/cc-pVDZ. Thus, only when dynamical correlation is included in a multi-reference method, the D 2h structure becomesthelowestinenergy. As follows from Table 2.1, using different optimized geometries has only a minor effect onΔE’s. For example, the energy differences calculated at the CCSD optimized geometries are within 0.03 eV from those calculated using the structures optimized at thesameleveloftheory,andtheobservedtrendsarethesame. 75 Thus, the presented ΔE’s exhibit the same trend as the optimized geometries and charge localization patterns — the more correlation we add the more stable the D 2h structurebecomes. Energy profiles from Fig. 2.8 give more detailed information on the shape of the PES. Due to the HF instabilities, the CCSD and CCSD(T) curves are discontinuous, unless the symmetry-broken UHF reference is used at D 2h . Apart from the cusp, the shapes of the OO-CCD and CCSD curves are surprisingly similar. The difference becomes more profound when triples are included: the OO-CCD(T) curve smoothly descendstowardstheD 2h minimum,whiletheCCSD(T)curveisgraduallyclimbingup all the way till the final drop at D 2h . Even though the corresponding D 2h → C 2v ΔE’s are close, the shape of the EOM-SF-CCSD and OO-CCD(T) curves are quite different: the former is much flatter at D 2h than the latter. The EOM-EA-CCSD curve parallels EOM-SF-CCSD around D 2h , but soon becomes spoiled by the instability in the triplet reference. Using optimized orbitals triplet reference solves the problem of HF instabil- ityandtheEOM-EA-ODcurveisverysimilartotheEOM-SF-CCSDone. TheB3LYP curveissteeperthantheEOM-SFone,consistentwithlargevalueofthecorresponding frequency(seenextsection). Overall,Fig. 2.8demonstratesthatdescribingtheshapeof PES is more challenging than calculating energy difference between two minima. One shouldexpectconsiderabledifferencesinthecorrespondingvibrationalfrequency. Note that the analytic CCSD or CCSD(T) frequencies would differ dramatically from those calculated by finite differences: the analytic frequencies will characterize the curvature ofthePEScorrespondingtothesymmetry-purereference,whereasfinitedifferencecal- culationswillsamplethecusp. To conclude, the analysis of the equilibrium structures, charge distributions, and energydifferencescalculatedatdifferentlevelsoftheorysuggeststhattheC 2v minimum is an artifact of incomplete treatment of electron correlation and HF instabilities in the 76 doublet or triplet references. The strong vibronic interactions result in rather flat PES along the D 2h → C 2v distortion. In the next section, we present our calculations of the photoelectron spectrum, which allows us to assess the quality of calculated geometries andtounambiguouslydeterminethesymmetryoftheanion. 2.5 Vibrationalfrequenciesandthephotoelectronspec- trumofpara-C 6 H − 4 Photoelectron spectroscopy is a very sensitive structural tool. The spectra depend on equilibriumgeometries,frequenciesandnormalmodesoftheinitial(radicalanion)and target(diradical)states. Theequilibriumgeometrieswerediscussedintheprevioussec- tion. The vibrational frequencies of the anion (at D 2h ) and the diradical are given in Table 2.2. The CCSD and B3LYP frequencies of the anion appear to be similar, except forthelowestb 1u modethatdescribesD 2h →C 2v displacementsandismostaffectedby vibronic interactions. This mode also exhibits the largest change relative to the neu- tral, i.e., it is considerably softer in C 6 H − 4 . The B3LYP frequency is 40% higher than the CCSD value, in favor of our assumption that B3LYP underestimates the magni- tudeofvibronicinteraction. Unfortunately,theonlysignificantvibrationalprogressions present the Franck-Condon factors (see section 2.3) are those for fully-symmetric nor- mal modes that exhibit displacements relative to the anion. Modes with significant fre- quencychange,suchastheaboveb 1u mode,yieldmuchsmallerfeaturesthatarelikely tobeburiedundermajorprogressions,asithappensforthepara-benzyneanion. The electron photodetachment spectrum of the para-benzyne radical anion reported by Wenthold and coworkers 12 consists of the two overlapping bands corresponding to thesinglet(X 1 A g )andtriplet( 3 B u )statesofthediradicalwiththeorigins,i.e.,electron 77 Table2.2: Vibrationalfrequenciesofthepara-benzyneradicalanionandthediradical. Symmetry X 2 A 1 a X 2 A 1 b a 3 B 1u a X 1 A 1 c a g 3137 3085 3221 3317 a g 1471 1434 1565 1440 a g 1189 1175 1166 1202 a g 992 981 1028 1065 a g 624 631 615 660 a u 879 917 920 994 a u 380 398 368 450 b 1g 691 721 810 722 b 2g 1120 890 871 961 b 2g 647 601 437 616 b 3g 3114 3061 3207 3302 b 3g 1579 1549 1652 1738 b 3g 1297 1284 1300 1327 b 3g 608 612 585 597 b 1u 3111 3056 3206 3300 b 1u 1445 1425 1471 1511 b 1u 1071 1060 1046 1101 b 1u 452 639 947 973 b 2u 3136 3082 3220 3316 b 2u 1331 1332 1375 1448 b 2u 1189 1203 1267 1254 b 2u 1027 1020 1081 1076 b 3u 715 718 753 796 b 3u 422 447 400 477 a CalculatedatUHF-CCSD/6-311+G** b CalculatedatB3LYP/6-311+G** c FromRef.39,SF-5050/6-31G** affinities (EAs), at 1.265± 0.008 eV and 1.430± 0.015 eV, respectively. Long vibra- tional progressions for frequencies at 635± 20cm −1 and 990± 20cm −1 (assigned to the singlet band), at 610± 15cm −1 and 995± 20cm −1 (assigned to the triplet band), 78 as well as a hotband at 615 ± 30 cm −1 were reported. For both the singlet and the tripletbands,thelower(about600cm −1 )andthehigher(about1000cm −1 )energypro- gressionswereassignedtosymmetricringdeformationandringbreathing,respectively. These bands, four lowest a g vibrations, are also present in all the computed spectra. Including the lowestb 1u mode resulted in very small features, completely obscured by thedenselinesofmajorprogressions. In the original paper 12 , the normal displacements were derived from the fitting pro- cedure. Table 2.3 and Fig. 2.9 compare these values with the shifts computed in this workfromthecorrespondingD 2h optimizedgeometriesandnormalmodesoftheanion. ThetripletandthesingletstatesofthediradicalaredescribedbytheCCSDandspin-flip DFTmethods,respectively(seesection2.3). ThecomputedCCSDshiftsareveryclose to those derived from the experiment. Values obtained with the CCSD normal modes and different equilibrium structures [e.g., CCSD(T), EOM-EA-CCSD, OO-CCD, and B3LYP] are quite similar to each other. However, the full B3LYP shifts (the last line in Table 2.3), i.e., those computed using the B3LYP normal modes and B3LYP equi- librium structures, are much larger for the ring breathing mode, which will result in a longerprogression. Finally,Fig. 2.10comparesthecalculatedandexperimentalphotodetachmentspec- tra. The upper panel of Fig. 2.10 shows the singlet and triplet bands of the spectrum for theD 2h structures calculated using the CCSD description of the anion. To compare with the experimental spectrum, the scaled intensities of the singlet and triplet bands wereaddedtogether. Twoscalingcoefficientsweredeterminedfromthefollowingcon- ditions: (i) the combined spectrum is normalized to 1; and (ii) the intensity of peak number 3 in the experimental and the calculated spectra coincide. The experimental spectrum was normalized to 1 as well. The bottom part of Fig. 2.10 compares the computed D 2h and the experimental spectra. The agreement between the two spectra 79 Table2.3: Displacementsinnormalmodesuponionization,A∗ √ amu. Method ν 5 a ν 4 b ν 5 a ν 4 b Singlet Triplet Exp. fit c ] -0.44 -0.14 0.76 -0.12 CCSD d -0.49 -0.21 -0.69 -0.09 CCSD(T) d -0.50 -0.27 -0.71 -0.14 EOM-EA-CCSD d -0.49 -0.22 -0.70 -0.10 OO-CCD d -0.49 -0.21 -0.69 -0.09 DFT/B3LYP d -0.47 -0.17 -0.67 -0.05 DFT/B3LYP e -0.40 -0.29 -0.60 -0.23 a Ringdeformation b Ringbreathing c Ref.12 d usingtheCCSDnormalmodes e usingtheDFT/B3LYPnormalmodes S -Q 6 3 5 S -Q 9 9 0 T -Q 6 1 0 T -Q 9 9 5 -0 .8 -0 .6 -0 .4 -0 .2 0 .0 e x p . fit C C S D C C S D (T ) E O M -E A -C C S D O O -C C D D F T D F T -fu ll D is p la c e m e n ts , A * a m u 1 /2 N o rm a l m o d e s Figure2.9: Displacementsalongthenormalmodesuponionization(D 2h structure). 80 1 .2 1 .4 1 .6 1 .8 2 .0 2 .2 2 .4 In te n s ity E le c tro n b in d in g e n e rg y , e V 1 .2 1 .4 1 .6 1 .8 2 .0 2 .2 2 .4 In te n s ity E le c tro n b in d in g e n e rg y , e V Figure2.10: Top: Singlet(solidline)andtriplet(dottedline)bandsoftheelectronpho- todetachmentspectrumforD 2h usingtheCCSDequilibriumgeometriesandfrequencies oftheanion. Bottom: theexperimental(solidline)andthecalculated(dottedline)spec- traforD 2h usingtheCCSDstructureandnormalmodes. 81 is excellent: the computed spectrum has the same features as the experimental one and differs only slightly in the peaks’ intensities. Fig. 2.11 shows the spectrum calculated for the C 2v CCSD structure, which is markedly different — it features additional pro- gressions that are not present in the experimental spectrum. They correspond to the ring deformation band and the combination bands involving this mode. To conclude, an excellent agreement between the spectrum calculated using the D 2h structure and the experimental one, along with significant differences between D 2h and C 2v , unam- biguouslyprovestheD 2h structureofthepara-benzyneradicalanion. Moreover,italso confirms highquality of the CCSDequilibrium structure andfrequencies for theanion, as well as the accuracy of the spin-flip description of the diradical. We also calculated the spectrum using the B3LYP optimized geometries and frequencies for p-C 6 H − 4 (and thesamestructuresandfrequenciesforthediradical,asabove). Theresultingspectrum (Fig. 2.12) has significantly different ratio of intensities for the two main bands and much longer progressions. Thus, even though DFT correctly predicts symmetry of the para-benzyneradicalanion,theoverallstructureisofpoorquality. Thereasonsforthis arediscussedbelow. 2.6 DFTself-interactionerrorandvibronicinteractions The analysis of equilibrium structures, charge localization, and energy differences between D 2h and C 2v structures presented in section 2.4 demonstrated that the magni- tude of symmetry breaking decreases when higher level methods are employed. More- over, the PES scans and the analysis of UHF references strongly suggest that the exis- tence of a C 2v minimum is an artifact. Indeed, EOM-SF-CCSD, the only CC method 82 1 .2 1 .4 1 .6 1 .8 2 .0 2 .2 2 .4 In te n s ity E le c tro n b in d in g e n e rg y , e V Figure 2.11: The experimental (solid line) and the calculated (dotted line) spectra for C 2v usingtheCCSDstructureandnormalmodes. thatemploys awell-behaved UHFreference (thatis, symmetry-pure, notstronglyspin- contaminated and stable), yields only a D 2h minimum and a smooth PES. The EOM- EA-OD potential energy profile parallels the EOM-SF one. Finally, modeling of the photoelectron spectrum ruled out the C 2v structure. Thus, we conclude that symmetry breakinginp-C 6 H − 4 ispurelyartefactual. DFT/B3LYP appears to be more robust with respect to this artefactual symmetry breaking as it yields only aD 2h structure. The important question is why, and whether ornotoneshouldconsiderDFTareliabletoolformodelingsystemswithstrongvibronic interactions. Detailedbenchmarkstudies(e.g.,Ref.18)providemanyexamplesofDFT giving “the right answer” for difficult symmetry breaking systems. Stanton explained theseobservationsbymakingconnectionsbetweenDFTresponsetheory,TD-DFTexci- tation energies and their Tamm-Dancoff approximation (footnote 23 in Ref. 20), and 83 1 .2 1 .4 1 .6 1 .8 2 .0 2 .2 2 .4 In te n s ity E le c tro n b in d in g e n e rg y , e V Figure 2.12: The experimental (solid line) and the calculated (dotted line) spectra for D 2h usingtheB3LYPstructureandnormalmodes. concluded that DFT describes pseudo Jahn-Teller effects reasonably well “for the right reason”, that is, because the poles appear in approximately the right places. Other examples, however, suggested that DFT might actually underestimate the magnitude ofvibronicinteractions 21 . Tocomplicatemattersevenfurther,wewouldliketomention the dehydro-meta-xylylene anion (DMX − ) for which DFT yields low symmetry (C s or C 2 ) equilibrium geometries 26 , while wave function based methods predict planar C 2v structures 45 . To gain insight, we begin by analyzing the character of interacting electronic states inp-C 6 H − 4 andDMX − atsymmetricandsymmetry-brokengeometries. Asdescribedin sections2.2and2.4,thetwointeractingstatesofp-C 6 H − 4 havetheexcesschargeequally delocalized between two radical centers (C 1 and C 4 ) at D 2h , whereas C 2v distortions 84 result in the charge localization at one of the carbons. DMX − shows exactly the oppo- sitebehavior 26,45 : attheplanarC 2v geometries,mostoftheexcesschargeresidesonthe σ-radical center, whereas C 2 or C s displacements facilitate the flow of the charge into the π-system. Thus, in both cases B3LYP favors the structures with more delocalized excess charge, which immediately brings about the infamous H + 2 example 46 . While HF dissociation curve for a one-electron H + 2 is exact, many DFT methods result in quasi- bound potential energy curves and with errors as large as 50 kcal/mol! This happens because DFT over-stabilizes the solutions with an electron being equally split between the two hydrogens, which is an artifact of self-interaction error (SIE) present in many functionals 47–49 . SIE is considerably larger in systems with odd number of electrons. Self-interaction energy, Coulomb interaction of an electron with itself, is always posi- tive,andelectrondelocalizationlowersitsmagnitudethuscausingartificialstabilization ofdelocalizedconfigurations. Otherwellunderstoodexamplesofartificialstabilization of delocalized states due to SIE include underestimated barriers for radical reactions, dissociationofcationicradicals,andmixed-valencetransitionmetaldimers 49 . As shown in Ref. 48, the magnitude of SIE in H + 2 strongly depends on the distance between the atoms, varying from a negligibly small value around the equilibrium to 55 kcal/mol (B3LYP) at the dissociation limit. The charge-bearing radical centers in the para-benzyne radical anion are 2.9 ˚ A apart. At this distance, the H + 2 curve exhibits as much as 60% of its maximum SIE. This suggests substantial SIE in the para-benzyne radicalanion. To test this assumption, we computed energy differences between the D 2h and C 2v minima (using the CCSD optimized geometries) employing functionals with a varying fraction of HF exchange. The results are shown in Fig. 2.13. The trend is clear: the relative stabilization of theD 2h structure is inversely proportional to the fraction of the HF exchange in the functional, and, consequently, is proportional to SIE. Combining 85 different correlation functionals with the fixed fraction of the HF exchange has little effect on the ΔE’s, for example, all functionals with 0% HF exchange yield approxi- matly the same D 2h - C 2v enegry separition. Therefore, we conclude that for DFT, the strengths of vibronic interactions and, consequently, a relative ordering of the D 2h and C 2v structures, are governed by SIE. Thus, B3LYP yields the right structure for the “wrong reason”. Therefore, it is not at all surprising that the quality of the B3LYP structure and the shape of PES is relatively poor, as follows from the comparison of the computed photo-electron spectrum with the experimental one. In DMX − , B3LYP yieldsthewrongstructure,butforthesame“rightreason”(SIE).Therefore,B3LYP(or any functional with SIE) is not a reliable method for the systems with strong vibronic interactions that affect charge localization patters, or when the interacting states are characterizedbyconsiderablydifferentchargedistributions. 0 2 0 4 0 6 0 8 0 1 0 0 -1 .0 -0 .8 -0 .6 -0 .4 -0 .2 0 .0 0 .2 H F V W N L Y P 7 5 -2 5 5 0 -5 0 B 3 L Y P B 3 L Y P 5 B 3 P W 9 1 E D F 2 E D F 1 G ill B L Y P S V W N B V W N G G 9 9 E (C 2 v ) - E (D 2 h ), e V H F e x c h a n g e , % Figure 2.13: Energy differences between the D 2h and C 2v minima (using the CCSD structures)calculatedusingdifferentdensityfunctionals. 86 2.7 Calculationofelectronaffinitiesofp-benzyne Different strategies for calculating electron affinities in problematic open-shell systems were discussed in Ref. 26. The simplest ”brute-force” approach, i.e., to take the differ- encebetweenthetotalenergiesoftheanionandtheneutralcalculatedatthesamelevel of theory, can give accurate values only if the errors in the corresponding total energies are very small (which can be achieved only at very high levels of theory), or if both species are described on the equal footing by the method and errors in the total ener- gies cancel out. The latter is difficult to achieve in the case when both the neutral and anionwavefunctionshaveacomplicatedopen-shellcharacter. Forexample,thesinglet stateofthediradicalissignificantlymulti-configurationalanditwouldbedescribedless accuratelythantheanionbysinglereferencemethods. Thetripletstateofthediradical, however, is single-configurational and can be described by single-reference methods with approximately the same accuracy as the doublet. Thus, triplet EAs calculated by energydifferencesemployingcorrelatedsingle-referencemethodsaremorereliablethan singlet EAs. The latter can be accurately computed by subtracting the experimental (or a calculated by an appropriate method) singlet—triplet gap from the calculated triplet electron affinity. This approach was used for calculating EAs of triradicals 26,45 and diradicals 50 . Belowwedemonstratethatthisschemeinthespiritofisodesmicreactions indeedyieldsaccuratesingletEAs. EAscalculatedatdifferentlevelsoftheoryforthesingletandthetripletstatesofthe diradical are presented in Table 2.4 and Fig. 2.14. The first two columns present triplet and singlet EAs computed by energy differences. The last column presents singlet EA calculated from the corresponding triplet EA and the best theoretical estimate of the singlet-triplet gap. For the triplet state, which is predominantly single-configurational, even relatively low level methods using a moderate basis, e.g. B3LYP and MP2 in 6- 311+G**,yieldreasonableEAs,e.g.,within0.11eVand0.07eVfromtheexperiment, 87 respectively. Forthesingletstate,however,EAscalculatedbyDFTandMP2are1.03eV and0.67eVoff,becausethesemethodsarenotappropriateforthemulti-configurational wave function of the singlet. Note that such brute-force approach reverses the states ordering in the diradical even at the CCSD level. For both states coupled-cluster meth- ods, especially with triples’ corrections, are in a reasonable agreement with the exper- imental values. For example, ROHF-CCSD(T) is within 0.16 eV for triplet and within 0.05eVforsinglet. Table 2.4: Electron affinities (eV) of the a 3 B 1u and X 1 A g states of the para-benzyne diradical. Method TripletEA SingletEA SingletEA a Experiment 1.430 1.265 Experiment-ZPE 1.375 1.140 6-311+G**basis B3LYP 1.481 2.172 1.310 UHF -0.994 2.679 -1.165 UHF-MP2 1.547 0.467 1.376 UHF-CCSD 0.893 1.659 0.722 UHF-CCSD(T) 1.207 1.085 1.036 ROHF -1.148 2.356 -1.319 ROHF-MP2 2.053 0.914 1.882 ROHF-CCSD 0.872 1.629 0.701 ROHF-CCSD(T) 1.216 1.094 1.045 OO-CCD 0.875 1.644 0.704 OO-CCD(T) 1.210 1.080 1.039 aug-cc-pVTZbasis ROHF -1.199 2.297 -1.370 ROHF-MP2 2.266 1.010 2.095 ROHF-CCSD 0.998 1.766 0.827 ROHF-CCSD(T) 1.395 1.258 1.224 a CalculatedbyusingtripletEAandthesinglet-tripletgapfromRef.51. 88 0 .5 1 .0 1 .5 2 .0 2 .5 6 -3 1 1 + G * * a u g -c c -p V T Z R O H F -C C S D R O H F -C C S D (T ) E A , e V U H F D F T /B 3 L Y P U H F -M P 2 U H F -C C S D (T ) U H F -C C S D O O -C C D O O -C C D (T ) E x p . -1 .0 -0 .5 0 .0 0 .5 1 .0 1 .5 2 .0 6 -3 1 1 + G * * a u g -c c -p V T Z R O H F -C C S D R O H F -C C S D (T ) E A , e V U H F D F T /B 3 L Y P U H F -M P 2 U H F -C C S D (T ) U H F -C C S D O O -C C D O O -C C D (T ) E x p . Figure2.14: Adiabaticelectronaffinitiescalculatedbyenergydifferencesforthesinglet (toppanel)andtriplet(bottompanel)stateofthepara-benzynediradical. 89 Forquantitativethermochemicalresults,largerbasissetsarerequired. Weperformed additionalcoupled-clustercalculationswiththeaug-cc-pVTZbasis(seeTable2.4). EA for the triplet calculated by CCSD(T) is within 0.02 eV from the experimental one, as expected for this method. However, the singlet CCSD(T) value is 0.12 eV off, because ofthepoordescriptionofthemulti-configurationalsingletstatebyground-statesingle- referencemethods. By employing the best theoretical estimate for the singlet-triplet gap in para- benzyne, i.e., 0.171 eV from Ref. 51, singlet EA calculated from the CCSD(T) triplet valueiswithin0.08eVfromtheexperiment,ascomparedto0.12eVforthebrute-force approach. AsfollowsfromtheexcellentagreementoftripletEA,theaccuracyofsinglet EAcanbeimprovedbyrefiningthevalueofthesinglet-tripletgap. 2.8 Conclusions Theshapeofground-statePESofC 6 H − 4 isstronglyaffectedbythevibronicinteractions with a low-lying excited state. The magnitude of vibronic interactions depends cru- cially on the the energy gap and couplings between the interacting states. Complicated open-shell character of the wave functions of the interacting states and HF instabilities challenge ab initio methodology: the shape of PES along the vibronic coupling coordi- natediffersdramaticallyfordifferentmethods. Forexample,manywavefunctionbased methodspredictvibronicinteractionstobestrongenoughtoresultinalower-symmetry C 2v minimum, in addition to the D 2h structure. EOM-SF-CCSD and B3LYP predict onlyasingleminimum,however,thecurvatureofthesurfaceisquitedifferent. The degree of this symmetry lowering can be quantified by structural parameters (e.g., the magnitude of nuclear displacements), charge localization and the resulting dipole moment, as well as energy difference between the two minima. We found that 90 the magnitude of symmetry breaking (as characterized by the above metrics) decreases when higher-level methods are employed. For example, the dipole moment of the C 2v structure decreases in the HF→CCSD→EOM-EA-CCSD series. Likewise, coupled- clustermethodswithtriplescorrectionsaswellasEOM-EA/SF-CCSDpredictthat D 2h isalowerenergystructure. Moreover,EOM-SFpredictedonly D 2h minima. Calculation of the potential energy profiles along the D 2h → C 2v distortion along with the stability analysisofUHFreferencesallowedustoattributethelower-symmetryminimumtothe incompletetreatmentofelectroncorrelationandUHFinstabilities. ThequalityofD 2h structurewasassessedbycomparingthecomputedphotoelectron spectrum against the experiment. The spectrum calculated using the CCSD description of the anion and the SF-DFT results for the diradical is in an excellent agreement with the experimental one, which is particularly encouraging in view of its dense and com- plicated nature. The spectrum calculated using the C 2v geometry and frequencies is markedlydifferent,whichallowsustoruleouttheC 2v structure. Aspointedoutinearlierstudies 16 ,B3LYPappearstoberobustw.r.t. thisartefactual symmetrybreaking,andyieldsonlyD 2h structure. Theanalysisofthestructureandthe potentialenergyprofilealongtheD 2h →C 2v distortionrevealsthatthestructureandthe curvature are considerably different from those calculated by reliable wave functions methods. Moreover, the computed spectrum differs significantly from the experimental one. Theshapeofthesurface,aswellasthecorrespondingharmonicfrequencyreveals that B3LYP underestimates the magnitude of vibronic interactions, as pointed out ear- lier by Crawford and coworkers 21 . Further analysis allowed us to attribute the failure of B3LYP to the self-interaction error. A similar behavior (that yielded seemingly dif- ferent result) was observed in DMX − , where SIE resulted in B3LYP overestimation of vibronic interactions and, consequently, a lower symmetry structure. In both cases, 91 large SIE originates in significant changes in charge localization patterns along sym- metry breaking coordinate. Thus, B3LYP is not a reliable method for systems where vibronicinteractionsarecoupledtodifferentchargelocalizationpatterns,ashappensin distonicradicalanddiradicalanions. TheworkpresentedinthischapterwaspublishedinRef.52. 92 Chapter2References [1] Wenk,H.H.;Winkler,M.;Sander,W.Angew.Chem.Int.Ed.Engl.2003,42,502. [2] Jones,R.R.;Bergman,R.G.J.Am.Chem.Soc.1972,94,660. [3] Lee, M.D. ; Dunne, T.S. ; Siegel, M.M. ; Chang, C.C. ; Morton, G.0. ; Borders, D.B.J.Am.Chem.Soc.1987,109,3464. [4] Borders,D.B.,Doyle,T.W.,Eds.Enediyneantibioticsasantitumoragents; Mar- celDeekker,NewYork,1995. [5] Maeda, H. , Edo, K. , Ishida, N. , Eds. Neocarzinostatin: the past, present, and futureofananticancerdrug; Springer,NewYork,1997. [6] Schottelius,M.J.;Chen,P.Angew.Chem.Int.Ed.Engl.1996,3,1478. [7] Hoffner, J.; Schottelius, M.J.; Feichtinger, D.; Chen, P. J. Am. Chem. Soc.1998, 120,376. [8] Nicolaou,K.C.;Smith,A.L.Acc.Chem.Res.1992,25,497. [9] Bowles,D.M.;Palmer,G.J.;Landis,C.A.;Scott,J.L.;Anthony,J.E.Tetrahedron 2001,57,3753. [10] Zeidan,T.;Manoharan,M.;Alabugin,I.V.J.Org.Chem.2006,71,954. [11] Wenthold,P.G.;Hu,J.;Squires,R.R.J.Am.Chem.Soc.1996,118,11865. [12] Wenthold, P.G. ; Squires, R.R. ; Lineberger, W.C. J. Am. Chem. Soc. 1998, 120, 5279. [13] Eshdat, L. ; Berger, H. ; Hopf, H. ; Rabinovitz, M. J. Am. Chem. Soc. 2002, 124, 3822. [14] Alabugin,I.V.;Kovalenko,S.V.J.Am.Chem.Soc.2002,124,9052. [15] Alabugin,I.V.;Manoharan,M.J.Am.Chem.Soc.2003,125,4495. [16] Nash,J.J.;Squires,R.R.J.Am.Chem.Soc.1996,118,11872. [17] Davidson,E.R.;Borden,W.T.J.Phys.Chem.1983,87,4783. [18] Cohen,R.D.;Sherrill,C.D.J.Chem.Phys.2001,114,8257. 93 [19] Crawford, T.D. ; Stanton, J.F. ; Allen, W.D. ; Schaefer III, H.F. J. Chem. Phys. 1997,107,10626. [20] Stanton,J.F.J.Chem.Phys.2001,115,10382. [21] Russ,N.J.;Crawford,T.D.;Tschumper,G.S.J.Chem.Phys.2004,120,7298. [22] L¨ owdin,P.-O. Rev.Mod.Phys.1963,35,496. [23] Eisfeld,W.;Morokuma,K.J.Chem.Phys.2000,113,5587. [24] Crawford, T.D. ; Kraka, E. ; Stanton, J.F. ; Cremer, D. J. Chem. Phys. 2001, 114, 10638. [25] Wladyslawski,M.;Nooijen,M.InACSSymposiumSeries,Vol.828,pages65–92, 2002. [26] Slipchenko,L.V.;Krylov,A.I.J.Phys.Chem.A2006,110,291. [27] Krishnan,R.;Binkley,J.S.;Seeger,R.;Pople,J.A.J.Chem.Phys.1980,72,650. [28] Frisch,M.J.;Pople,J.A.;Binkley,J.S.J.Chem.Phys.1984,80,3265. [29] Woon,D.E.;Jr.,T.H.DunningJ.Chem.Phys.1994,100. [30] Cristian,A.M.C.;Shao,Y.;Krylov,A.I.J.Phys.Chem.A2004,108,6581. [31] Purvis,G.D.;Bartlett,R.J.J.Chem.Phys1982,76,1910. [32] Raghavachari, K. ; Trucks, G.W. ; Pople, J.A. ; Head-Gordon, M. Chem. Phys. Lett.1989,157,479. [33] Sherrill,C.D.;Krylov,A.I.;Byrd,E.F.C.;Head-Gordon,M. J.Chem.Phys.1998, 109,4171. [34] Nooijen,M.;Bartlett,R.J.J.Chem.Phys.1995,102,3629. [35] Krylov,A.I.Chem.Phys.Lett.2001,338,375. [36] Levchenko,S.V.;Krylov,A.I.J.Chem.Phys.2004,120,175. [37] Becke,A.D.J.Chem.Phys.1993,98,5648. [38] Arnold,D.W.Ph.D.Thesis,UCBerkeley,1994. [39] Shao,Y.;Head-Gordon,M.;Krylov,A.I. J.Chem.Phys.2003,118,4807. 94 [40] NBO 5.0. Glendening, E.D. ; Badenhoop, J.K. ; Reed, A.E. ; Carpenter, J.E. ; Bohmann, J.A. ; Morales, C.M. ; Weinhold, F. Theoretical Chemistry Institute, UniversityofWisconsin,Madison,WI,2001. [41] Mulliken,R.S.J.Chem.Phys.1955,23,1997. [42] Shao, Y. ; Molnar, L.F. ; Jung, Y. ; Kussmann, J. ; Ochsenfeld, C. ; Brown, S. ; Gilbert, A.T.B. ; Slipchenko, L.V. ; Levchenko, S.V. ; O’Neil, D. P. ; Distasio Jr., R.A.;Lochan,R.C.;Wang,T.;Beran,G.J.O.;Besley,N.A.;Herbert,J.M.;Lin, C.Y. ; Van Voorhis, T. ; Chien, S.H. ; Sodt, A. ; Steele, R.P. ; Rassolov, V. A. ; Maslen,P.;Korambath,P.P.;Adamson,R.D.;Austin,B.;Baker,J.;Bird,E.F.C. ; Daschel, H. ; Doerksen, R.J. ; Drew, A. ; Dunietz, B.D. ; Dutoi, A.D. ; Furlani, T.R. ; Gwaltney, S.R. ; Heyden, A. ; Hirata, S. ; Hsu, C.-P. ; Kedziora, G.S. ; Khalliulin, R.Z. ; Klunziger, P. ; Lee, A.M. ; Liang, W.Z. ; Lotan, I. ; Nair, N. ; Peters, B. ; Proynov, E.I. ; Pieniazek, P.A. ; Rhee, Y.M. ; Ritchie, J. ; Rosta, E. ; Sherrill,C.D.;Simmonett,A.C.;Subotnik,J.E.;WoodcockIII,H.L.;Zhang,W. ;Bell,A.T.;Chakraborty,A.K.;Chipman,D.M.;Keil,F.J.;Warshel,A.;Herhe, W.J.;SchaeferIII,H.F.;Kong,J.;Krylov,A.I.;Gill,P.M.W.;Head-Gordon,M. Phys.Chem.Chem.Phys.2006,8,3172. [43] ACES II. Stanton, J.F. ; Gauss, J. ; Watts, J.D. ; Lauderdale, W.J. ; Bartlett, R.J. 1993. The package also contains modified versions of the MOLECULE Gaussian integral program of J. Alml¨ of and P.R. Taylor, the ABACUS integral derivative program written by T.U. Helgaker, H.J.Aa. Jensen, P. Jørgensen and P.R. Taylor, andthePROPSpropertyevaluationintegralcodeofP.R.Taylor. [44] Helgaker,T.;Jørgensen,P.;Olsen,J.Molecularelectronicstructuretheory;Wiley &Sons,2000. [45] Munsch, T.E. ; Slipchenko, L.V. ; Krylov, A.I. ; Wenthold, P.G. J. Org. Chem. 2004,69,5735. [46] Bally,T.;Sastry,G.N.J.Phys.Chem.A1997,101,7923. [47] Polo,V.;Kraka,E.;Cremer,D.MolecularPhysics2002,100,1771. [48] Zhang,Y.;Yang,W.J.Chem.Phys.1998,109,2604. [49] Lundber,M.;Siegbahn,P.E.M.J.Chem.Phys.2005,122,1. [50] Reed,D.R.;Hare,M.;Kass,S.R.J.Am.Chem.Soc.2000,122,10689. [51] Slipchenko,L.V.;Krylov,A.I.J.Chem.Phys.2002,117,4694. [52] Vanovschi,V.;Krylov,A.I.;Wenthold,P.G.Theor.Chim.Acta2008,120,45. 95 Chapter3 Futurework Theeffectivefragmentpotentialmethod,whichemergedinthelastdecade,iscurrently anactiveareaofresearchinboththemethodologyandapplications. Theimplementation developed in this work increases the availability of the method and opens the possibil- ity to combine it with various ab initio theories available in Q-Chem 1 (e.g., RI-CCSD, EOM-SF-CCSD).Thus,itislikelytoincreasetheinteresttothismethodanditspoten- tial even further. Since EFP is a relatively new method, the opportunities for the future developments are abound. Below the most interesting directions of the future work are presentedintheseparatesectionsforthemethod,implementationandapplications. 3.1 Methodology Three-bodydispersion A recent work by Podeszwa ar al. 2 demonstrated that three-body dispersion can con- tribute up to 50% to total three-body interaction energy for non-polar molecular clus- ters, such as cyclic benzene timer. Therefore, we expect that three body dispersion will become increasingly important in the three-dimensional aromatic clusters and in the bulk. To model intermolecular interactions in such systems accurately, the three-body dispersion term needs to be included in the EFP model. The approach, which must not compromise the computational efficiency of EFP, is yet to be developed. One of the 96 possible directions is to apply the approximations made in the derivation of the two- body EFP dispersion 3 to the three-body dispersion expression in terms of frequency- dependentelectricpolarizabilitiesderivedinRef.4. Impovedelectrostaticsexpansion In the present EFP formulation 5 , the exchange-repulsion, polarization and dispersion pointsaretheLMOcentroids,whereastheelectrostaticsexpansionpointsarethenuclei and bond midpoints. The LMOs centroids account for both bond midpoints and lone pairs. Thus,thelonepairs,whichalsorepresentconcentrationoftheelectronicdensity, are not included in the electrostatic expansion. An interesting question is if using the LMO centroids and the nuclei (instead of the bond midpoints and the nuclei) as the electrostatics expansion points will make the EFP description more balanced and thus moreaccurate. Additionally,themultipolesattheLMOcentroidsmightbetterfitinthe unifiedEFPdampingschemeproposedinRef.29. LinearscalingEFP ThecomputationalcomplexityofallEFPcomponentsintheirstraightforwardformula- tion is O(N 2 ), where N is the number of fragments. Application of the fast multipole method 6 totheelectrostaticsandpolarizationcomponentsalongwithacutoffprocedure to fast decaying exchange-repulsion and dispersion, will bring the computational com- plexity down toO(Nlog(N)). This will open the possibility for modeling much larger molecularsystemsandthepropertiesofthebulkcondensedphaseattheEFPlevel. Correlatedtreatmentoftheactiveregion In the present EFP formulation 5 the active region is treated at the HF level of theory, which does not account for electron correlation. However, it is well established 7–18 that 97 accurate description of many molecular properties requires correlated methods. Elec- troncorrelationisespeciallyimportantfortheopen-shellspecies,suchasmulticonfigu- rationalexcitedstates,whichrequireinclusionofatleaststaticcorrelation. The implementation of EFP developed in this work is a part of Q-Chem. This state oftheart ab initioelectronicstructurepackageimplementsawidevarietyofcorrelated methods 1 . Therefore,apromisingdirectionofthefutureworkistodevelopacorrelated EFP theory and to extend our implementation so that the active region can be treated with any of correlation methods available in Q-Chem. This theory will allow one to account for the environmental effects in the electronic structure studies of radicals and electronically excited species, as well as to model chemical reactions in the condensed phase. Relevant to the research conducted in this work (chapter 2), correlated EFP will allowtomodelthereactionofBergmancyclizationinsolutions. Below a possible framework for such correlated EFP theory is outlined. In the presentEFPformulation,theeffectivefragmentsaffecttheactiveregionviaone-electron terms in the Hamiltonian. These terms originate from the electrostatics and polariza- tion components. Correlated methods can similarly use this modified Hamiltonian to account for the field of the effective fragments. However, a complication arises from the polarization component. Its contribution to the Hamiltonian depends on the values of the induced dipoles. The induced dipoles depend on the wavefunction of the active region, which in turn depends on the Hamiltonian. Thus the correlated wavefunction andtheinduceddipolesshouldbeconsistentwitheachother. Tosatisfytheconsistency requirement,theinduceddipolesandthecorrelatedwavefunctionsshouldbecomputed self-consistentlywiththefollowingiterativeprocedure: 1. Computethecorrelatedwavefunctionforagivenvaluesoftheinduceddipoles. 2. Updatetheinduceddipolestoaccountforthenewwavefunction. 98 3. If the induced dipoles and the wavefunction are converged exit, otherwise return tostep1. However, such strict self-consistency requires to recompute the correlated wavefunc- tionmultipletimesand,therefore,iscomputationallyexpensive. Alternatively,onemay use the induced dipoles consistent with the HF wavefunction, which is computed as a prerequisite for all single reference correlated methods. This way the consistency is compromisedtoreducethecomputationalcost. Yetanotherapproachistoperformjust one step toward the consistency, i.e., to compute the correlated wavefunction based on theHFinduceddipoles,updatetheinduceddipolesandrecomputethecorrelatedwave- function again. The accuracy of the last two approximations is yet to be investigated. For the methods computing the correlated wavefunction iteratively (such as CCSD), it ispossibletodesignamoreefficientintegratedapproach,similartothatfortheHartree- FockSCFprocedure(section1.3). Inthisapproach,ateveryiterationboththecorrelated wavefunction and the induced dipoles are updated and thus consistence of both will be achieved once they converge. The drawback is the increased complexity in the imple- mentation. 3.2 Implementation Gradients As discussed in section 1.7, the energy minima for the benzene dimer structures were obtainedusingscansofthepotentialenergysurface. Forlargesystems,abetterapproach would be to use a gradient-based optimization procedure. The gradient of the electro- static energy was formulated in the original EFP paper 19 . Gradient for the dispersion energyissimplyaderivativeof 1 R 6 . Moreintricateanddifficulttoimplementarerecently 99 publishedgradientsforpolarization 20,21 andexchange-repulsionenergies 22 . Implement- ing the gradient of EFP energy will allow one to perform geometry optimizations and moleculardynamicssimulationsforlargemolecularsystems. Fragmentparamters In the developed implementation, the EFP fragment parameters are defined in the Q- Chem standard input, as described in section 1.6. The parameters of an effective frag- ment are computed from an individual isolated molecules, as discussed in sections 1.2, 1.31.4and1.5. Thus,foragivenmoleculethesamesetofEFPparameterscanbeused in all computations. This immediately prompts creating a library of the pre-computed EFP parameters for molecules that frequently occur in the systems studied with EFP. Such library can, for instance, include the parameters for common organic solvents (e.g., water, alcohols, acetone, DMSO, acetonitrile, formic acid, dichloromethane, car- bontetrachloride,chloroform,tetrahydrofuran,benzene,phenol,toluene,etc.). TheEFP parameters included in the library can be carefully selected and benchmarked against high level ab initio theories. Therefore, the benefits of such a library are not just the convenienceitprovides,butalsoahigh-qualityassurancefortheEFPparameters. QM-EFdispersionandexchange-repulsioninteractions InthecurrentstateofEFPtheory,theinteractionbetweentheeffectivefragmentsandthe abinitioregionisduetotheelectrostaticsandpolarizationcomponentsonly. Thus,there is no account for dispersion and exchange-repulsion interactions between them. The proper theories for such interactions are presently developed in M.S. Gordon’s group. Meanwhile,itispossibletoaccountforsuchinteractionsbytreatingtheabinitiosystem as an effective fragment for these two types of interactions. Since the wavefunction of the ab initiopartisreadilyavailable,thedispersionandexchange-repulsionparameters 100 can be computed as discussed in sections 1.4 and 1.5, respectively. From these param- eters, dispersion and exchange-repulsion interaction energies between the active region andtheeffectivefragmentscanbeevaluated. Implementingsuchschemewillallowone to properly compute the binding energy in the systems with an ab initio region. In the present EFP theory, binding energy is properly computed for systems composed of the effective fragments only, whereas molecular systems with an ab initio region are used to model the solvation effects. Implementing the suggested technique will allow one to combineaccuratedescriptionofthesolvationeffectsandbindingenergies. 3.3 Applications Numerous applications of EFP for modeling the solvation effects and the binding ener- giesintheclustersandcondensedphasecanbeenvisioned. Here,thefuturedirectionof the research of multiple π−π interactions conducted in this work (section 1.7) is pre- sented. As demonstrated in Ref. 2 using the cyclic benzene trimer example, three-body dispersion has a substantial contribution to many-body interaction energy. Introduction of the three-body dispersion term in EFP is proposed in section 3.1. Once this term is available,theextendedEFPmethodcanbeappliedto: • studytheroleofmany-body π−π interactionsinthemultidimensionalaromatic clustersandthebulk; • model the mechanical properties of carbon-based materials, such as graphite and fullerenepowders; • investigatetheroleofπ−π stackingintherigidityofDNA. 101 Chapter3References [1] Shao, Y. ; Molnar, L.F. ; Jung, Y. ; Kussmann, J. ; Ochsenfeld, C. ; Brown, S. ; Gilbert, A.T.B. ; Slipchenko, L.V. ; Levchenko, S.V. ; O’Neil, D. P. ; Distasio Jr., R.A.;Lochan,R.C.;Wang,T.;Beran,G.J.O.;Besley,N.A.;Herbert,J.M.;Lin, C.Y. ; Van Voorhis, T. ; Chien, S.H. ; Sodt, A. ; Steele, R.P. ; Rassolov, V. A. ; Maslen,P.;Korambath,P.P.;Adamson,R.D.;Austin,B.;Baker,J.;Bird,E.F.C. ; Daschel, H. ; Doerksen, R.J. ; Drew, A. ; Dunietz, B.D. ; Dutoi, A.D. ; Furlani, T.R. ; Gwaltney, S.R. ; Heyden, A. ; Hirata, S. ; Hsu, C.-P. ; Kedziora, G.S. ; Khalliulin, R.Z. ; Klunziger, P. ; Lee, A.M. ; Liang, W.Z. ; Lotan, I. ; Nair, N. ; Peters, B. ; Proynov, E.I. ; Pieniazek, P.A. ; Rhee, Y.M. ; Ritchie, J. ; Rosta, E. ; Sherrill,C.D.;Simmonett,A.C.;Subotnik,J.E.;WoodcockIII,H.L.;Zhang,W. ;Bell,A.T.;Chakraborty,A.K.;Chipman,D.M.;Keil,F.J.;Warshel,A.;Herhe, W.J.;SchaeferIII,H.F.;Kong,J.;Krylov,A.I.;Gill,P.M.W.;Head-Gordon,M. Phys.Chem.Chem.Phys.2006,8,3172. [2] Podeszwa,R.;Szalewicz,K.J.Chem.Phys.2007,126,194101. [3] Adamovix,I.;Gordon,M.S.Mol.Phys.2005,103,379. [4] McLachlan,A.D.Mol.Phys.1963,6,423. [5] Gordon, M.S. ; Slipchenko, L. ; Li, H. ; Jensen, J.H. In Annual Reports in Com- putational Chemistry; Spellmeyer, D.C. , Wheeler, R. , Eds., Vol. 3 of Annual ReportsinComputationalChemistry; Elsevier,2007; pages177–193. [6] Schmidt,K.E.;Lee,M.A.J.Stat.Phys.1991,63,1223. [7] Woon,D.E.;Dunning,T.H.J.Chem.Phys.1993,99,1914. [8] Peterson,K.A.;Kendall,R.A.;Dunning,T.H.J.Chem.Phys.1993,99,1930. [9] Peterson,K.A.;Kendall,R.A.;Dunning,T.H.J.Chem.Phys.1993,99,9790. [10] Peterson,K.A.;Woon,D.E.;Dunning,T.H.J.Chem.Phys.1994,100,7410. [11] Woon,D.E.J.Chem.Phys.1994,100,2838. [12] Woon,D.E.;Dunning,T.H.J.Chem.Phys.1994,101,8877. [13] Peterson,K.A.;Dunning,T.H.J.Chem.Phys.1995,102,2032. [14] Peterson,K.A.;Dunning,T.H.J.Chem.Phys.1997,106,4119. [15] Woon,D.E.;Peterson,K.A.;Dunning,T.H.J.Chem.Phys.1998,109,2233. 102 [16] Wilson,A.K.;Dunning,T.H.J.Chem.Phys.1997,106,8718. [17] Peterson,K.A.;Dunning,T.H.J.Phys.Chem.A1997,101,6280. [18] Peterson, K.A.; Wilson, A.K.; Woon, D.E.; Dunning, T.H. Theor. Chem. Acc. 1997,97,251. [19] Day, P.N. ; Jensen, J.H. ; Gordon, M.S. ; Webb, S.P. ; Stevens, W.J. ; M.Krauss; D.Garmer;H.Basch;D.CohenJ.Chem.Phys.1996,105,1968. [20] Li,H.;Gordon,M.S.J.Chem.Phys.2007,126. [21] Li,H.;Netzloff,H.M.;Gordon,M.S.J.Chem.Phys.2006,125. [22] Li,H.;Gordon,M.S.Theor.Chem.Acc.2006,115,385. 103 Bibliography [1] I.AdamovicandM.S.Gordon.SolventeffectsontheS(N)2reaction: Application of the density functional theory-based effective fragment potential method. J. Phys.Chem.A,109(8):1629–1636,2005. [2] I.AdamovixandM.S.Gordon. Dynamicpolarizability,dispersioncoefficientC6 and dispersion energy in the effective fragment potential method. Mol. Phys., 103(2-3):379–387,2005. [3] I.V. Alabugin and S.V. Kovalenko. C1-C5 photochemical cyclization of enediynes. J.Am.Chem.Soc.,124:9052–9053,2002. [4] I.V. Alabugin and M. Manoharan. Radical-anionic cyclizations of enediynes: Remarkableeffectsofbenzannelationandremotesubstituentsoncycliorearoma- tizationreactions. J.Am.Chem.Soc.,125:4495–4509,2003. [5] N.L. Allunger, Y.H. Yuh, and J.H. Lii. Molecular mechanics - the MM3 force- fieldforhydrocarbons.1. J.Am.Chem.Soc.,111(23):8551–8566,1989. [6] D.W.Arnold. StudyofRadicals,Clusters,andTransitionStateSpeciesbyAnion PhotoelectronSpectroscopy. PhDthesis,UCBerkeley,1994. [7] B.Askew,P.Ballester,C.Buhr,K.S.Jeong,S.Jones,K.Parris,K.Williams,and J.Rebek. Molecularrecognitionwithconvergentfunctionalgroups.VI.Synthetic andstructuralstudieswithamodelreceptorfornucleicacidcomponents. J. Am. Chem.Soc.,111(3):1082–1090,1989. [8] R.Balawender,B.Safi,andP.Geerlings. Solventeffectontheglobalandatomic DFT-based reactivity descriptors using the effective fragment potential model. Solvationofammonia. J.Phys.Chem.A,105(27):6703–6710,2001. [9] T.BallyandG.N.Sastry.Incorrectdissociationbehaviorofradicalionsindensity functionalcalculations. J.Phys.Chem.A,101:7923–7925,1997. [10] A.D.Becke. Density-functionalthermochemistry.iii.Theroleofexactexchange. J.Chem.Phys.,98:5648,1993. 104 [11] D.B.BordersandT.W.Doyle,editors. Enediyneantibioticsasantitumoragents. MarcelDeekker,NewYork,1995. [12] D.R. Bowler, T. Miyazaki, and M.J. Gillan. Recent progress in linear scaling ab initio electronic structure techniques. J. Phys.: Condens. Matter, 14(11):2781– 2798,2002. [13] D.M. Bowles, G.J. Palmer, C.A. Landis, J.L. Scott, and J.E. Anthony. The Bergman reaction as a synthetic tool: Advantages and restrictions. Tetrahedron, 57:3753–3760,2001. [14] S.F.Boys. Constructionofsomemolecularorbitalstobeapproximatelyinvariant forchangesfromonemoleculetoanother.Rev.Mod.Phys.,32(2):296–299,1960. [15] S.F.Boys. InQuantumScienceofAtoms,Molecules,andSolids,pages253–262. AcademicPress,1966. [16] M.F. Brana, M. Cacho, A. Gradillas, B. Pascual-Teresa, and A. Ramos. Interca- latorsasanticancerdrugs. Curr.Pharm.Des.,7(17):1745–1780,2001. [17] A.D. Buckingham. Molecular quadrupole moments. Q. Rev. Chem. Soc., 13(3):183–214,1959. [18] S.K. Burley and G.A. Petsko. Aromatic-aromatic interaction: A mechanism of proteinstructurestabilization. Science,229(4708):23–28,1985. [19] C.J. Casewit, K.S. Colwell, and A.K. Rappe. Application of a universal force- fieldtoorganic-molecules. J.Am.Chem.Soc.,114(25):10035–10046,1992. [20] K.R.S. Chandrakumar, T.K. Ghanty, S.K. Ghosh, and T. Mukherjee. Hydra- tion of uranyl cations: Effective fragment potential approach. J. Molec. Struct. (Theochem),807(1-3):93–99,2007. [21] R.D. Cohen and C.D. Sherrill. The performance of density functional theory for equilibrium molecular properties of symmetry breaking molecules. J. Chem. Phys.,114:8257–8269,2001. [22] T.D. Crawford, E. Kraka, J.F. Stanton, and D. Cremer. Problematic p-benzyne: Orbital instabilities, biradical character, and broken symmetry. J. Chem. Phys., 114:10638–10650,2001. [23] T.D. Crawford, J.F. Stanton, W.D. Allen, and H.F. Schaefer III. Hartree-Fock orbitalinstabilityenvelopesinhighlycorrelatedsingle-referencewavefunctions. J.Chem.Phys.,107:10626–10632,1997. 105 [24] A.M.C. Cristian, Y. Shao, and A.I. Krylov. Bonding patterns in benzene trirad- icals from structural, spectroscopic, and thermochemical perspectives. J. Phys. Chem.A,108:6581–6588,2004. [25] E.R. Davidson and W.T. Borden. Symmetry breaking in polyatomic molecules: Realandartifactual. J.Phys.Chem.,87:4783–4790,1983. [26] P.N. Day, J.H. Jensen, M.S. Gordon, S.P. Webb, W.J. Stevens, M.Krauss, D.Garmer, H.Basch, and D.Cohen. An effective fragment method for modeling solventeffectsinquantummechanicalcalculations.J.Chem.Phys.,105(5):1968– 1986,1996. [27] P.N. Day, M.S. Pachter, R.and Gordon, and G.N. Merrill. A study of water clus- tersusingtheeffectivefragmentpotentialandMonteCarlosimulatedannealing. J.Chem.Phys.,112(5):2063–2073,2000. [28] P.N. Day and R. Pachter. A study of aqueous glutamic acid using the effective fragmentpotentialmethod. J.Chem.Phys.,107(8):2990–2999,1997. [29] G.R. Desiraju and A. Gavezzotti. From molecular to crystal structure; polynu- cleararomatichydrocarbons. JCSChem.Comm.,pages621–623,1989. [30] W. Eisfeld and K. Morokuma. A detailed study on the symmetry breaking an its effectonthepotentialsurfaceofNO 3 . J.Chem.Phys.,113:5587–5597,2000. [31] L. Eshdat, H. Berger, H. Hopf, and M. Rabinovitz. Anionic cyclization of a cross-conjugatedenediyne. J.Am.Chem.Soc.,124:3822–3823,2002. [32] P.M. Felker, P.M. Maxton, and M.W. Schaeffer. Nonlinear raman studies of weakly bound complexes and clusters in molecular beams. Chem. Rev., 94(7):1787–1805,1994. [33] T.K. Firman and C.R. Landis. Valence bond concepts applied to the molecular mechanicsdescriptionofmolecularshapes.4.Transitionmetalswithπ-bonds. J. Am.Chem.Soc.,123(47):11728–11742,2001. [34] M.A. Freitag and M.S. Gordon. Using the effective fragment-potential model to study the solvation of acetic and formic acids. Abstracts of Papers of the AmericanChemicalSociety,218:U511–U511,1999. 1. [35] M.A. Freitag, M.S. Gordon, J.H. Jensen, and W.J. Stevens. Evaluation of charge penetration between distributed multipolar expansions. J. Chem. Phys., 112(17):7300–7306,2000. 106 [36] M.J. Frisch, J.A. Pople, and J.S. Binkley. Self-consistent molecular orbital methods 25. Supplementary functions for Gaussian basis sets. J. Chem. Phys., 80:3265–3269,1984. [37] J.L. Gao and D.G. Truhlar. Quantum mechanical methods for enzyme kinetics. Annu.Rev.Phys.Chem.,53:467–505,2002. [38] E.Gazit.Self-assembledpeptidenanostructures: thedesignofmolecularbuilding blocks and their technological utilization. Chem. Soc. Rev., 36(8):1263–1269, 2007. [39] E.D. Glendening, J.K. Badenhoop, A.E. Reed, J.E. Carpenter, J.A. Bohmann, C.M.Morales,andF.Weinhold. NBO5.0. TheoreticalChemistryInstitute,Uni- versityofWisconsin,Madison,WI,2001. [40] V.Gogonea,L.M.Westerhoff,andJr.K.M.Merz. Quantummechanical/quantum mechanicalmethods.I.Adivideandconquerstrategyforsolvingtheschrodinger equation for large molecular systems using a composite density functional– semiempiricalHamiltonian. J.Chem.Phys.,113(14):5604–5613,2000. [41] M.S. Gordon, M.A. Freitag, P. Bandyopadhyay, J.H. Jensen, V. Kairys, and W.J. Stevens. Theeffectivefragmentpotentialmethod: AQM-basedMMapproachto modelingenvironmentaleffectsinchemistry. J.Phys.Chem.A,105(2):293–307, 2001. [42] M.S. Gordon, L. Slipchenko, H. Li, and J.H. Jensen. Chapter 10 the effective fragment potential: A general method for predicting intermolecular interactions. In D.C. Spellmeyer and R. Wheeler, editors, Annual Reports in Computational Chemistry, volume 3 of Annual Reports in Computational Chemistry, pages 177 –193.Elsevier,2007. [43] T. Helgaker, P. Jørgensen, and J. Olsen. Molecular electronic structure theory. Wiley&Sons,2000. [44] J.R. Hill and J. Sauer. Molecular mechanics potential for silica and zeolite catalysts based on ab initio calculations. 2. Aluminosilicates. J. Phys. Chem., 99(23):9536–9550,1995. [45] J. Hoffner, M.J. Schottelius, D. Feichtinger, and P. Chen. J. Am. Chem. Soc., 120:376,1998. [46] C.A.HunterandJ.K.M.Sanders. Thenatureofπ−π interactions. J.Am.Chem. Soc.,112(14):5525–5534,1990. [47] ISO/IEC.ISO/IEC14882:2003: Programminglanguages: C++.Technicalreport, InternationalOrganizationforStandardization,Geneva,Switzerland.,2003. 107 [48] J.H. Jensen and M.S. Gordon. An approximate formula for the intermolecular Pauli repulsion between closed shell molecules. Mol. Phys., 89(5):1313–1325, 1996. [49] J.H. Jensen and M.S. Gordon. An approximate formula for the intermolecular Pauli repulsion between closed shell molecules. II. Application to the effective fragmentpotentialmethod. J.Chem.Phys.,108(12):4772–4782,1998. [50] R.R. Jones and R.G. Bergman. p-Benzyne. Generation as an intermediate in a thermal isomerization reaction and trapping evidence for the 1,4-benzenediyl structure. J.Am.Chem.Soc.,94:660,1972. [51] W.L. Jorgensen and N.A. McDonald. Development of an all-atom force field for heterocycles. Properties of liquid pyridine and diazenes. J. Molec. Struct. (Theochem),424(1-2):145–155,1998. [52] K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M. Uebayasi. Fragment molecu- lar orbital method: An approximate computational method for large molecules. Chem.Phys.Lett.,313(3-4):701–706,1999. [53] R. Krishnan, J.S. Binkley, R. Seeger, and J.A. Pople. Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions. J. Chem. Phys., 72:650,1980. [54] A.I.Krylov. Size-consistentwavefunctionsforbond-breaking: Theequation-of- motionspin-flipmodel. Chem.Phys.Lett.,338:375–384,2001. [55] M.D.Lee,T.S.Dunne,M.M.Siegel,C.C.Chang,G.0.Morton,andD.B.Borders. Calichemicins, a novel family of antitumor antibiotics. 1. Chemistry and partial structureofcalichemicin. J.Am.Chem.Soc.,109:3464,1987. [56] K.R. Leopold, G.T. Fraser, S.E. Novick, and W. Klemperer. Current themes in microwave and infrared spectroscopy of weakly bound complexes. Chem. Rev., 94(7):1807–1827,1994. [57] S.V. Levchenko and A.I. Krylov. Equation-of-motion spin-flip coupled-cluster model with single and double substitutions: Theory and application to cyclobu- tadiene. J.Chem.Phys.,120(1):175–185,2004. [58] H. Li and M.S. Gordon. Gradients of the exchange-repulsion energy in the gen- eral effective fragment potential method. Theor. Chem. Acc., 115(5):385–390, 2006. [59] H. Li and M.S. Gordon. Polarization energy gradients in combined quantum mechanics,effectivefragmentpotential,andpolarizablecontinuummodelcalcu- lations. J.Chem.Phys.,126(12),2007. 108 [60] H. Li, M.S. Gordon, and J.H. Jensen. Charge transfer interaction in the effective fragmentpotentialmethod. J.Chem.Phys.,124(21),2006. [61] H. Li, H.M. Netzloff, and M.S. Gordon. Gradients of the polarization energy in theeffectivefragmentpotentialmethod. J.Chem.Phys.,125(19),2006. [62] H. Lin and D.G. Truhlar. QM/MM: what have we learned, where are we, and wheredowegofromhere? Theor.Chem.Acc.,117(2):185–199,2007. [63] G.Lippert,J.Hutter,andM.Parrinello.Ahybridgaussianandplanewavedensity functionalscheme. Mol.Phys.,92(3):477–487,1997. [64] P.-O.L ¨ owdin. Rev.Mod.Phys.,35:496,1963. [65] M. Lundber and P.E.M. Siegbahn. Quantifying the effects of the self-interaction error in DFT: When do the delocalized states appear? J. Chem. Phys., 122(224103):1–9,2005. [66] A.D.MacKerell,D.Bashford,Bellott,R.L.Dunbrack,J.D.Evanseck,M.J.Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F.T.K. Lau, C. Mattos, S. Michnick, T. Ngo, D.T. Nguyen, B. Prodhom, W.E. Reiher, B. Roux, M. Schlenkrich, J.C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz-Kuczera, D. Yin, and M. Karplus. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B, 102(18):3586–3616,1998. [67] H. Maeda, K. Edo, and N. Ishida, editors. Neocarzinostatin: the past, present, andfutureofananticancerdrug. Springer,NewYork,1997. [68] S.L.Mayo,B.D.Olafson,andW.A.Goddard. DREIDING-agenericforce-field formolecularsimulations. J.Phys.Chem.,94(26):8897–8909,1990. [69] A.D.McLachlan. Three-bodydispersionforces. Mol.Phys.,6(4):423–427,1963. [70] G.N. Merrill and M.S. Gordon. Study of small water clusters using the effective fragmentpotentialmodel. J.Phys.Chem.A,102(16):2650–2657,1998. [71] G.N. Merrill and S.P. Webb. Anion-water clusters A(-)(H2O)(1-6), A = OH, F, SH, Cl, and Br. an effective fragment potential test case. J. Phys. Chem. A, 107(39):7852–7860,2003. [72] G.N. Merrill and S.P. Webb. The application of the effective fragment potential method to molecular anion solvation: A study of ten oxyanion-water clusters, A(-)(H2O)(1-4). J.Phys.Chem.A,108(5):833–839,2004. 109 [73] G.N.Merrill,S.P.Webb,andD.B.Bivin. Formationofalkalimetal/alkalineearth cation water clusters, M(H2O)(1-6), M = Li+, Na+, K+, Mg2+, and Ca2+: An effectivefragmentpotential(EFP)casestudy. J.Phys.Chem.A,107(3):386–396, 2003. [74] G. Monard and K.M. Merz. Combined quantum mechanical/molecular mechan- ical methodologies applied to biomolecular systems. Acc. Chem. Res., 32(10):904–911,1999. [75] K. Mueller-Dethlefs, O. Dopfer, and T.G. Wright. ZEKE spectroscopy of com- plexesandclusters. Chem.Rev.,94(7):1845–1871,1994. [76] K. Muller-Dethlefs and P. Hobza. Noncovalent interactions: A challenge for experimentandtheory. Chem.Rev.,100(1):143–168,2000. [77] R.S. Mulliken. Electronic population analysis on LCAO-MO molecular wave functions.I. J.Chem.Phys.,23(10):1833–1840,1955. [78] R.S. Mulliken. Report on notation for the spectra of polyatomic molecules. J. Chem.Phys.,23:1997,1955. [79] T.E. Munsch, L.V. Slipchenko, A.I. Krylov, and P.G. Wenthold. Reactivity and structure of the 5-dehydro-m-xylylene anion. J. Org. Chem., 69:5735–5741, 2004. [80] J.J.NashandR.R.Squires. Theoreticalstudiesofo-,m-,andp-benzynenegative ions. J.Am.Chem.Soc.,118:11872–11883,1996. [81] H.M.NetzloffandM.S.Gordon. Theeffectivefragmentpotential: Smallclusters andradialdistributionfunctions. J.Chem.Phys.,121(6):2711–2714,2004. [82] K.C. Nicolaou and A.L. Smith. Molecular design, chemical synthesis, and bio- logicalactionofenediynes. Acc.Chem.Res.,25:497–503,1992. [83] M. Nooijen and R.J. Bartlett. Equation of motion coupled cluster method for electronattachment. J.Chem.Phys.,102:3629–3647,1995. [84] K. A. Peterson and T. H. Dunning. Benchmark calculations with correlated molecular wave-functions. 7. Binding-energy and structure of the hf dimer. J. Chem.Phys.,102(5):2032–2041,1995. [85] K. A. Peterson and T. H. Dunning. Benchmark calculations with correlated molecular wave functions. 11. Energetics of the elementary reactions F+H-2, O+H-2,andH’+HCl. J.Phys.Chem.A,101(35):6280–6292,1997. 110 [86] K. A. Peterson and T. H. Dunning. Benchmark calculations with correlated molecular wave functions. 8. Bond energies and equilibrium geometries of the chnandc2hn(n=1-4)series. J.Chem.Phys.,106(10):4119–4140,1997. [87] K.A.Peterson, R.A.Kendall, andT.H.Dunning. Benchmarkcalculationswith correlatedmolecularwave-functions.2.Configuration-interactioncalculationson 1st-rowdiatomichydrides. J.Chem.Phys.,99(3):1930–1944,1993. [88] K.A.Peterson, R.A.Kendall, andT.H.Dunning. Benchmarkcalculationswith correlatedmolecularwave-functions.3.Configuration-interactioncalculationson 1strowhomonucleardiatomics. J.Chem.Phys.,99(12):9790–9805,1993. [89] K.A.Peterson,A.K.Wilson,D.E.Woon,andT.H.Dunning. Benchmarkcalcu- lations with correlated molecular wave functions. 12. Core correlation effects on the homonuclear diatomic molecules b-2-f-2. Theor. Chem. Acc., 97(1-4):251– 259,1997. [90] K. A. Peterson, D. E. Woon, and T. H. Dunning. Benchmark calculations with correlatedmolecularwave-functions.4.Theclassicalbarrierheightoftheh+h-2- ]h-2+hreaction. J.Chem.Phys.,100(10):7410–7415,1994. [91] R. Podeszwa. Comment on ”Beyond the benzene dimer: An investigation of the additivityofπ−π interactions”. J.Phys.Chem.A,112(37):8884–8885,2008. [92] R. Podeszwa, R. Bukowski, and K. Szalewicz. Potential energy surface for the benzene dimer and perturbational analysis ofπ−π interactions. J. Phys. Chem. A,110(34):10345–10354,2006. [93] R. Podeszwa and K. Szalewicz. Three-body symmetry-adapted perturbation theory based on kohn-sham description of the monomers. J. Chem. Phys., 126(19):194101,2007. [94] V. Polo, E. Kraka, and D. Cremer. Electron correlation and the self-interaction errorofdensityfunctionaltheory. MolecularPhysics,100:1771–1790,2002. [95] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press,NewYork,NY,USA,2007. [96] G.D. Purvis and R.J. Bartlett. A full coupled-cluster singles and doubles model: Theinclusionofdisconnectedtriples. J.Chem.Phys,76:1910,1982. [97] J.R. Rabinowitz, K. Namboodiri, and H. Weinstein. A finite expansion method for the calculation and interpretation of molecular electrostatic potentials. Int. J. Quant.Chem.,29(6):1697–1704,1986. 111 [98] K. Raghavachari, G.W. Trucks, J.A. Pople, and M. Head-Gordon. A fifth-order perturbation comparison of electron correlation theories. Chem. Phys. Lett., 157:479–483,1989. [99] D.R.Reed,M.Hare,andS.R.Kass.Formationofgas-phasedianionsanddistonic ions as a general method for the synthesis of protected reactive intermediates. Energetics of 2,3- and 2,6-dehydronaphthalene. J. Am. Chem. Soc., 122:10689– 10696,2000. [100] R.C.RizzoandW.L.Jorgensen. OPLSall-atommodelforamines: Resolutionof theaminehydrationproblem. J.Am.Chem.Soc.,121(20):4827–4836,1999. [101] G.V. Rossum. The Python Language Reference Manual. Network Theory Ltd., 2003. [102] N.J.Russ,T.D.Crawford,andG.S.Tschumper. Realversusartifactualsymmetry breaking effects in Hartree-Fock, density-functional, and coupled-cluster meth- ods. J.Chem.Phys.,120:7298–7306,2004. [103] W. Saenger. Principles of Nucleic Acid Structure. Springer-Verlag: New York, 1984. [104] A.Sanz-GarciaandR.J.Wheatley. Atomicrepresentationofthedispersioninter- actionenergy. ChemPhysChem,5(5):801–807,2003. [105] K.E. Schmidt and M.A. Lee. Implementing the fast multipole method in three dimensions. J.Stat.Phys.,63(5):1223–1235,1991. [106] M.W. Schmidt, K.K. Baldridge, S.T. Elbert J.A. Boatz, M.S. Gordon, S. Koseki J.H. Jensen, N. Mastunaga, K.A. Nguyen, T.L. Windus S. Su, M. Dupuis, and J.A. Montgomery. General atomic and molecular electronic structure system. J. Comput.Chem.,14(11):1347–1363,1993. [107] M.J.SchotteliusandP.Chen. Angew.Chem.Int.Ed.Engl.,3:1478,1996. [108] M.Sch¨ utz,G.Hetzer,andH.Werner.Low-orderscalinglocalelectroncorrelation methods.I.LinearscalinglocalMP2. J.Chem.Phys.,111(13):5691–5705,1999. [109] E. Schwegler and M. Challacombe. Linear scaling computation of the Hartree- Fockexchangematrix. J.Chem.Phys.,105(7):2726–2734,1996. [110] G.E. Scuseria and P.Y. Ayala. Linear scaling coupled cluster and perturbation theoriesintheatomicorbitalbasis. J.Chem.Phys.,111(18):8330–8343,1999. [111] Y. Shao, M. Head-Gordon, and A.I. Krylov. The spin-flip approach within time- dependent density functional theory: Theory and applications to diradicals. J. Chem.Phys.,118:4807–4818,2003. 112 [112] Y. Shao, L.F. Molnar, Y. Jung, J. Kussmann, C. Ochsenfeld, S. Brown, A.T.B. Gilbert, L.V. Slipchenko, S.V. Levchenko, D. P. O’Neil, R.A. Distasio Jr., R.C. Lochan, T. Wang, G.J.O. Beran, N.A. Besley, J.M. Herbert, C.Y. Lin, T. Van Voorhis,S.H.Chien,A.Sodt,R.P.Steele,V.A.Rassolov,P.Maslen,P.P.Koram- bath, R.D. Adamson, B. Austin, J. Baker, E.F.C. Bird, H. Daschel, R.J. Doerk- sen, A. Drew, B.D. Dunietz, A.D. Dutoi, T.R. Furlani, S.R. Gwaltney, A. Hey- den, S. Hirata, C.-P. Hsu, G.S. Kedziora, R.Z. Khalliulin, P. Klunziger, A.M. Lee, W.Z. Liang, I. Lotan, N. Nair, B. Peters, E.I. Proynov, P.A. Pieniazek, Y.M. Rhee, J. Ritchie, E. Rosta, C.D. Sherrill, A.C. Simmonett, J.E. Subotnik, H.L. WoodcockIII,W.Zhang,A.T.Bell,A.K.Chakraborty,D.M.Chipman,F.J.Keil, A.Warshel,W.J.Herhe,H.F.SchaeferIII,J.Kong,A.I.Krylov,P.M.W.Gill,and M. Head-Gordon. Advances in methods and algorithms in a modern quantum chemistryprogrampackage. Phys.Chem.Chem.Phys.,8:3172–3191,2006. [113] C.D.Sherrill,A.I.Krylov,E.F.C.Byrd,andM.Head-Gordon. Energiesandana- lytic gradients for a coupled-cluster doubles model using variational Brueckner orbitals: Application to symmetry breaking in O + 4 . J. Chem. Phys., 109:4171, 1998. [114] M.SierkaandJ.Sauer. Findingtransitionstructuresinextendedsystems: Astrat- egy based on a combined quantum mechanics-empirical valence bond approach. J.Chem.Phys.,112(16):6983–6996,2000. [115] M.O.SinnokrotandC.D.Sherrill. High-accuracyquantummechanicalstudiesof π−π interactions in benzene dimers. J. Phys. Chem. A, 110(37):10656–10668, 2006. [116] M.O.Sinnokrot,E.F.Valeev,andC.D.Sherrill. Estimatesoftheabinitiolimitfor π−πinteractions: Thebenzenedimer. J.Am.Chem.Soc.,124(36):10887–10893, 2002. [117] L.V.SlipchenkoandM.S.Gordon. Electrostaticenergyintheeffectivefragment potential method: Theory and application to benzene dimer. J. Comput. Chem., 28(1):276–291,2007. [118] L.V. Slipchenko and M.S. Gordon. Damping functions in the effective fragment potentialmethod. Mol.Phys.,2009. [119] L.V.SlipchenkoandA.I.Krylov. Singlet-tripletgapsindiradicalsbythespin-flip approach: Abenchmarkstudy. J.Chem.Phys.,117:4694–4708,2002. [120] L.V. Slipchenko and A.I. Krylov. Efficient strategies for accurate calculations of electronic excitation and ionization energies: Theory and application to the dehydro-meta-xylyleneanion. J.Phys.Chem.A,110:291–298,2006. 113 [121] T. Smith, L.V. Slipchenko, and M.S. Gordon. Modelingπ−π interactions with the effective fragment potential method: The benzene dimer and substituents. J. Phys.Chem.A,112(23):5286–5294,2008. [122] D.B. Smithrud and F. Diederich. Strength of molecular complexation of apolar solutes in water and in organic solvents is predictable by linear free energy rela- tionships: Ageneralmodelforsolvationeffectsonapolarbinding. J.Am.Chem. Soc.,112(1):339–343,1990. [123] J.F.Stanton.Coupled-clustertheory,pseudo-Jahn-Tellereffectsandconicalinter- sections. J.Chem.Phys.,115:10382–10393,2001. [124] J.F. Stanton, J. Gauss, J.D. Watts, W.J. Lauderdale, and R.J. Bartlett. ACES II, 1993. ThepackagealsocontainsmodifiedversionsoftheMOLECULEGaussian integral program of J. Alml¨ of and P.R. Taylor, the ABACUS integral derivative program written by T.U. Helgaker, H.J.Aa. Jensen, P. Jørgensen and P.R. Taylor, andthePROPSpropertyevaluationintegralcodeofP.R.Taylor. [125] A.J.Stone. Distributedmultipoleanalysis,orhowtodescribeamolecularcharge distribution. Chem.Phys.Lett.,83(2):233–239,1981. [126] A.J. Stone. In The Theory of Intermolecular Forces, pages 50 – 62. Oxford Uni- versityPress,1997. [127] A.J.Stone. Distributedmultipoleanalysis: Stabilityforlargebasissets. J.Chem. TheoryComput.,1(6):1128–1132,2005. [128] A.J.StoneandM.Alderton. Distributedmultipoleanalysis-methodsandappli- cations. Mol.Phys.,56(5):1047–1064,1985. [129] K.T. Tang and J.P. Toennies. An improved simple model for the van der Waals potentialbasedonuniversaldampingfunctionsforthedispersioncoefficients. J. Chem.Phys.,80(8):3726–3741,1984. [130] T.P. Tauer and C.D. Sherrill. Beyond the benzene dimer: An investigation of the additivityofπ−π interactions. J.Phys.Chem.A,109(46):10475–10478,2005. [131] S. Tsuzuki, T. Uchimaru, K. Matsumura, M. Mikami, and K. Tanabe. Effects of the higher electron correlation correction on the calculated intermolecular inter- action energies of benzene and naphthalene dimers: comparison between MP2 andCCSD(T)calculations. Chem.Phys.Lett.,319(5-6):547–554,2000. [132] V. Vanovschi, A.I. Krylov, and P.G. Wenthold. Structure, vibrational frequen- cies,ionizationenergies,andphotoelectronspectrumofthepara-benzyneradical anion. Theor.Chim.Acta,120:45–58,2008. 114 [133] H.H.Wenk,M.Winkler,andW.Sander. Onecenturyofarynechemistry. Angew. Chem.Int.Ed.Engl.,42:502–528,2003. [134] P.G. Wenthold, J. Hu, and R.R. Squires. J. Am. Chem. Soc., 118:11865–11871, 1996. [135] P.G. Wenthold, R.R. Squires, and W.C. Lineberger. Ultraviolet photoelectron spectroscopy of the o-, m-, and p-benzyne negative ions. Electron affinities and singlet-tripletsplittingsforo-,m-,andp-benzyne. J.Am.Chem.Soc.,120:5279– 5290,1998. [136] C.A. White, B.G. Johnson, P.M.W. Gill, and M. Head-Gordon. Linear scaling density functional calculations via the continuous fast multipole method. Chem. Phys.Lett.,253(3-4):268–278,1996. [137] A.K.WilsonandT.H.Dunning. Benchmarkcalculationswithcorrelatedmolec- ular wave functions. 10. Comparison with ”exact” mp2 calculations on ne, hf, h2o,andn-2. J.Chem.Phys.,106(21):8718–8726,1997. [138] M.WladyslawskiandM.Nooijen. ThephotoelectronspectrumoftheNO 3 radi- cal revisited: A theoretical investigation of potential energy surfaces and conical intersections. InACSSymposiumSeries,volume828,pages65–92,2002. [139] M.WolfsbergandL.Helmholz. Thespectraandelectronicstructureofthetetra- hedralionsMnO4-,CrO4–,andClO4-. J.Chem.Phys.,20(5):837–843,1952. [140] D. E. Woon. Benchmark calculations with correlated molecular wave-functions. 5. The determination of accurate ab-initio intermolecular potentials for he-2, ne- 2,andar-2. J.Chem.Phys.,100(4):2838–2850,1994. [141] D. E. Woon and T. H. Dunning. Benchmark calculations with correlated molec- ular wave-functions. 1. Multireference configuration-interaction calculations for the2nd-rowdiatomichydrides. J.Chem.Phys.,99(3):1914–1929,1993. [142] D. E. Woon and T. H. Dunning. Benchmark calculations with correlated molecularwave-functions.6.2nd-row-A(2)andfirst-row2nd-row-ABdiatomic- molecules. J.Chem.Phys.,101(10):8877–8893,1994. [143] D. E. Woon and T. H. Dunning Jr. Gaussian basis sets for use in correlated molecular calculations. iv. calculation of static electrical response properties. J. Chem.Phys.,100,1994. [144] D. E. Woon, K. A. Peterson, and T. H. Dunning. Benchmark calculations with correlated molecular wave functions. IX. The weakly bound complexes Ar-H-2 andAr-HCl. J.Chem.Phys.,109(6):2233–2241,1998. 115 [145] X.Ye,Z.Li,W.Wang,K.Fan,W.Xu,andZ.Hua. Theparallelπ−π stacking: amodelstudywithMP2andDFTmethods. Chem.Phys.Lett.,397(1-3):56–61, 2004. [146] S. Yoo, F. Zahariev, S. Sok, and M.S. Gordon. Solvent effects on opticalproper- tiesofmolecules: Acombinedtime-dependentdensityfunctionaltheory/effective fragmentpotentialapproach. J.Chem.Phys.,129(14):8,2008. [147] T. Zeidan, M. Manoharan, and I.V. Alabugin. Ortho effect in the Bergman cyclization: Interception of p-benzyne intermediate by intramolecular hydrogen abstraction. J.Org.Chem.,71:954–961,2006. [148] Y.ZhangandW.Yang. Achallengefordensityfunctionals: Self-interactionerror increases for systems with a noninteger number of electrons. J. Chem. Phys., 109(7):2604–2608,1998. 116 AppendixA GAMESSinputforgeneratingEFP parametersofbenzene TheGAMESSinputfileforgenerationtheEFPparametersattheHF/6-311G++(3df,2p) leveloftheoryforthebenzenemonomerusedinthisworkisgivenbelow. $CONTRL RUNTYP=MAKEFP ISPHER=1 $END $SYSTEM TIMLIM=99999 MWORDS=150 $END $BASIS GBASIS=N311 NGAUSS=6 NFFUNC=1 NDFUNC=3 NPFUNC=2 DIFFSP=.T. DIFFS=.T. $END $SCF SOSCF=.F. DIIS=.T. $END $STONE BIGEXP=4.0 $END $DAMP IFTTYP(1)=3,2 IFTFIX(1)=1,1 THRSH=500 $END $DAMPGS C3=C1 H4=H2 C5=C1 H6=H2 C7=C1 H8=H2 C9=C1 H10=H2 C11=C1 H12=H2 BO43=BO21 BO53=BO31 BO65=BO21 BO75=BO31 BO87=BO21 BO97=BO31 BO109=BO21 BO111=BO31 BO119=BO31 BO1211=BO21 $END $DATA benzene 117 C1 C1 6.0 -1.3942300 0.00000000 0.0000000 H2 1.0 -2.4764870 0.00000000 0.0000000 C3 6.0 -0.6971150 0.00000000 1.2074388 H4 1.0 -1.2382435 0.00000000 2.1447002 C5 6.0 0.6971150 0.00000000 1.2074388 H6 1.0 1.2382435 0.00000000 2.1447002 C7 6.0 1.3942300 0.00000000 -0.0000000 H8 1.0 2.4764870 0.00000000 0.0000000 C9 6.0 0.6971150 0.00000000 -1.2074388 H10 1.0 1.2382435 0.00000000 -2.1447002 C11 6.0 -0.6971150 0.00000000 -1.2074388 H12 1.0 -1.2382435 0.00000000 -2.1447002 $END 118 AppendixB Geometriesofbenezeneoligomers Geometries of all the benzeneoligomers studied in this work are provided below inthe format of the Q-Chem $efp fragments section. Fig. 1.16 provides the definition of the format. Sandwich(S)structureofbenzenedimer: benzene 0.0 0.0 0.0 benzene 0.0 4.0 0.0 T-shaped(T)structureofbenzenedimer: benzene 0.0 0.0 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 5.2 0.0 1.57079633 1.57079633 0.00000000 Parallel-displaced(PD)structureofbenzenedimer: benzene 0.0 0.0 0.0 benzene 1.2 3.9 0.0 Sbenzenetrimer: benzene 0.0 0.0 0.0 benzene 0.0 4.0 0.0 benzene 0.0 8.0 0.0 T1benzenetrimer: benzene 0.0 0.0 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 5.2 0.0 1.57079633 1.57079633 0.00000000 benzene 0.0 10.4 0.0 1.57079633 1.57079633 1.57079633 T2benzenetrimer: benzene 0.0 0.0 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 5.2 0.0 1.57079633 1.57079633 0.00000000 benzene 0.0 -5.2 0.0 1.57079633 1.57079633 0.00000000 PDbenzenetrimer: 119 benzene 0.0 0.0 0.0 benzene 1.2 3.9 0.0 benzene 0.0 7.8 0.0 Cbenzenetrimer: benzene 2.6 0.0 0.0 -0.25943951 0.00000000 0.00000000 benzene -2.6 0.0 0.0 1.83495559 0.00000000 0.00000000 benzene 0.0 -4.503332 0.0 3.92935069 0.00000000 0.00000000 PD-Tbenzenetrimer: benzene 0.0 0.0 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 5.2 0.0 1.57079633 1.57079633 0.00000000 benzene 0.0 -3.9 1.2 1.57079633 1.57079633 1.57079633 PD-Sbenzenetrimer: benzene 0.0 -4.0 0.0 benzene 0.0 0.0 0.0 benzene 1.2 3.9 0.0 S-Tbenzenetrimer: benzene 0.0 0.0 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 5.2 0.0 1.57079633 1.57079633 0.00000000 benzene 0.0 -4.0 0.0 1.57079633 1.57079633 1.57079633 S-10benzeneoligomer(largerS-oligomersareconstructedsimilarly): benzene 0.0 0.0 0.0 benzene 0.0 4.0 0.0 benzene 0.0 8.0 0.0 benzene 0.0 12.0 0.0 benzene 0.0 16.0 0.0 benzene 0.0 20.0 0.0 benzene 0.0 24.0 0.0 benzene 0.0 28.0 0.0 benzene 0.0 32.0 0.0 benzene 0.0 36.0 0.0 T-10benzeneoligomer(largerT-oligomersareconstructedsimilarly): benzene 0.0 0.0 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 5.2 0.0 1.57079633 1.57079633 0.00000000 benzene 0.0 10.4 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 15.6 0.0 1.57079633 1.57079633 0.00000000 benzene 0.0 20.8 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 26.0 0.0 1.57079633 1.57079633 0.00000000 benzene 0.0 31.2 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 36.4 0.0 1.57079633 1.57079633 0.00000000 benzene 0.0 41.6 0.0 1.57079633 1.57079633 1.57079633 benzene 0.0 46.8 0.0 1.57079633 1.57079633 0.00000000 120 PD-10benzeneoligomer(largerPD-oligomersareconstructedsimilarly): benzene 0.0 0.0 0.0 benzene 1.2 3.9 0.0 benzene 0.0 7.8 0.0 benzene 1.2 11.7 0.0 benzene 0.0 15.6 0.0 benzene 1.2 19.5 0.0 benzene 0.0 23.4 0.0 benzene 1.2 27.3 0.0 benzene 0.0 31.2 0.0 benzene 1.2 35.1 0.0 121 AppendixC ConverteroftheEFPparameters Ascript,whichconvertstheGAMESSinputparameterstotheQ-Chemformatisgiven below: #!/usr/bin/env python """EFP parameters converter. Converts from GAMESS to Q-Chem input format. Usage: python converter.py <GAMESS INPUT> """ import sys import re # Bohrs to Angstroms convertion factor B2A = 0.529177249 __version__ = "1.15" __author__ = "Vitalii Vanovschi" if "set" not in globals(): from sets import Set as set (FRAG, COORDS, MONOPOLES, DIPOLES, QUADRUPOLES, OCTUPOLES, \ BASIS, WAVEFUNCTION, FOCK, LMO, MULTIPLICITY, SCREEN2,\ POLARIZABLE, DYNAMIC, CANONVEC, CANONFOK, SCREEN3) = xrange(0, 17) def dump_qchem_frag(frag_name, labels, charges, dipoles, quadrupoles, \ octupoles, wavefunction, fock, centroids, screen2, \ ordered_labels, nuc_charges, pols, disps, basis, coordinates): """Outputs EFP fragment parameters in the Q-Chem format""" print "fragment", frag_name # output atoms for coordinate in coordinates: if coordinate[0][0:2] != "bo": print coordinate[0].rstrip("0123456789") + " " \ + " ".join(coordinate[1:]) # output atomic charges (as multipoles) for label in ordered_labels: coord = labels[label] if label in nuc_charges and float(nuc_charges[label]) > 0.1: print "mult ", " ".join(coord) print nuc_charges[label] # output multipoles for label in ordered_labels: 122 coord = labels[label] print "mult ", " ".join(coord) if label in screen2: print "cdamp " + screen2[label][1] if label in charges: print charges[label] if label in dipoles: print " ".join(dipoles[label]) if label in quadrupoles: print " ".join(quadrupoles[label]) if label in octupoles: print " ".join(octupoles[label]) # output polarization points for pol in pols: print "pol ", pol[1], " ", pol[2], " ", pol[3] print " ".join(pol[4:]) # output dispersion points dispindexes = list(set([x[0] for x in disps])) dispindexes.sort() for i in dispindexes: isfirst = True for disp in disps: if disp[0] == i: if isfirst: print "disp ", disp[1], " ", disp[2], " ", disp[3] isfirst = False print disp[4], print # output exchange-repulsion print "er_basis" print "er_wavefunction" wfindexes = wavefunction.keys() wfindexes.sort() if len(wfindexes) * (len(wfindexes) + 1) != 2 * len(fock): raise AssertionError("Inconsistent input.\n" "The following assertion has failed:\n" "len(wfindexes) * (len(wfindexes) + 1) == 2 * len(fock)\n" "Values: len(wfindexes) = " + str(len(wfindexes)) + ", " "len(fock) = " + str(len(fock)) + ".\n" "Did you remove ’>’ from the last line of the FOCK MATRIX \n" "section in GAMESS input? (an old GAMESS bug)") for ind in wfindexes: orb = wavefunction[ind] i = 0 for shell in basis: if shell == ’s’: i += 1 if shell == ’p’: i += 3 if shell == ’sp’ or shell == ’l’: i += 4 if shell == ’d’: tmp = orb[i:i+6] #10 -> 10 #11 -> 12 #12 -> 15 #13 -> 11 #14 -> 13 123 #15 -> 14 orb[i] = tmp[0] orb[i+1] = tmp[3] orb[i+2] = tmp[1] orb[i+3] = tmp[4] orb[i+4] = tmp[5] orb[i+5] = tmp[2] i += 6 if shell == ’f’: tmp = orb[i:i+10] orb[i] = tmp[0] #fxxx orb[i+1] = tmp[3] #fxxy orb[i+2] = tmp[5] #fxyy orb[i+3] = tmp[1] #fyyy orb[i+4] = tmp[4] #fxxz orb[i+5] = tmp[9] #fxyz orb[i+6] = tmp[6] #fyyz orb[i+7] = tmp[7] #fxzz orb[i+8] = tmp[8] #fyzz orb[i+9] = tmp[2] #fzzz i += 10 for i in range(len(orb)): if i % 5 == 0: print ind, if orb[i][0] == "-": print orb[i], else: print " " + orb[i], if i % 5 == 4: print print print "er_fock_matrix" for i in range(len(fock)): print fock[i], if i % 5 == 4: print print print "er_lmos" for lmo in centroids: print " ".join(lmo) def parse_gamess(gamess_file): """Parses GAMESS input file""" str2int = {"frag" : FRAG, "coordinates" : COORDS, "monopoles" : MONOPOLES, "dipoles" : DIPOLES, "quadrupoles" : QUADRUPOLES, "octupoles" : OCTUPOLES, "basis" : BASIS, "wavefunction" : WAVEFUNCTION, "fock" : FOCK, "lmo" : LMO, "multiplicity" : MULTIPLICITY, "screen2": SCREEN2, "polarizable": POLARIZABLE, "dynamic": DYNAMIC, "canonvec": CANONVEC, "canonfok": CANONFOK, 124 "screen3": SCREEN3 } frag_name = None ordered_labels = [] labels = {} coordinates = [] charges = {} nuc_charges = {} dipoles = {} quadrupoles = {} octupoles = {} screen2 = {} centroids = [] wavefunction = {} fock = [] expect = FRAG prefix = "" pols = [] disps = [] basis = [] print "$efp_params" for line1 in gamess_file: line_orig = line1 line1 = line1.strip() if not line1 or line1.lower() == ’STOP’: continue if line1[-1] == ">": prefix += " " + line1 continue else: line = prefix + " " + line1 prefix = "" fragment_start = re.search("\$([A-Za-z0-9]+)", line) if fragment_start: text = fragment_start.groups()[0].lower() if text in ("efrag", "contrl"): continue if text == "end" and frag_name: # print "frag_name:", frag_name # print "labels:", labels # print "dipoles:", dipoles # print "quadrupoles:", quadrupoles # print "octupoles:", octupoles # print "wavefunction:", wavefunction # print "fock:", fock # print "centroids:", centroids dump_qchem_frag(frag_name, labels, charges, dipoles, quadrupoles, octupoles, wavefunction, fock, centroids, screen2, ordered_labels, nuc_charges, pols, disps, basis, coordinates) frag_name = None labels = {} charges = {} nuc_charges = {} dipoles = {} 125 quadrupoles = {} octupoles = {} screen2 = {} centroids = [] coordinates = [] ordered_labels = [] wavefunction = {} fock = [] expect = FRAG prefix = "" pols = [] disps = [] basis = [] if text != "end": frag_name = text while 1: line = gamess_file.next() words = re.findall("([A-Za-z0-9\.+-]+)", line) if words[0].lower() == "coordinates": break expect = COORDS continue words = re.findall("([A-Za-z0-9\.+-]+)", line) if words: firstword = words[0].lower() else: firstword = None if firstword in str2int: expect = str2int[firstword] continue if firstword == ’projection’ and len(words) > 1: expect = str2int[words[1].lower()] continue if expect == COORDS and len(words) > 3: atom_x = str(float(words[1]) * B2A) atom_y = str(float(words[2]) * B2A) atom_z = str(float(words[3]) * B2A) labels[firstword] = (atom_x, atom_y, atom_z) coordinates.append((firstword, atom_x, atom_y, atom_z)) ordered_labels.append(firstword) if expect == MONOPOLES and len(words) > 1: charges[firstword] = words[1] nuc_charges[firstword] = words[2] if expect == DIPOLES and len(words) == 4: dipoles[firstword] = words[1:] if expect == QUADRUPOLES and len(words) == 7: quadrupoles[firstword] = words[1:] if expect == OCTUPOLES and len(words) == 11: octupoles[firstword] = words[1:] if expect == WAVEFUNCTION: # 3 1-1.40747045E-01 2.10323654E-01 2.42989533E-01-4.68907985E-01-1.03447318E-01 key = int(line_orig[:3].strip()) 126 if key not in wavefunction: wavefunction[key] = [] values = line_orig[5:].rstrip() while values: wavefunction[key].append(values[:15].strip()) values = values[15:] if expect == FOCK: fock += words if expect == LMO and len(words) == 4: centroids.append([str(float(x) * B2A) for x in words[1:]]) if expect == SCREEN2: screen2[firstword] = words[1:] if expect == POLARIZABLE: if len(words) == 4: pol_x = str(float(words[1]) * B2A) pol_y = str(float(words[2]) * B2A) pol_z = str(float(words[3]) * B2A) pols.append([words[0], pol_x, pol_y, pol_z]) elif len(words) == 9: pols[-1].extend(words) if expect == DYNAMIC: if len(words) == 5 or len(words) > 5 and words[5] == "--": index = int(words[1]) disp_x = str(float(words[2]) * B2A) disp_y = str(float(words[3]) * B2A) disp_z = str(float(words[4]) * B2A) disps.append([index, disp_x, disp_y, disp_z]) elif len(words) == 9: disps[-1].append(sum([float(x) for x in words[:3]]) / 3) if expect == BASIS: if firstword in (’s’, ’sp’, ’p’, ’d’, ’l’, ’f’): basis.append(firstword) print "$end" if __name__ == "__main__": if len(sys.argv) != 2: print "Usage: python converter.py <GAMESS INPUT>" sys.exit() INPUT = open(sys.argv[1]) parse_gamess(INPUT) INPUT.close() 127 AppendixD Fragment-fragmentelectrostatic interaction The equations for electrostatics inter-fragment interaction energy are presented below. k and l are the multipole expansion points. q is the charge, μ is the dipole, Θ is the quadrupole, Ω is the octupole, R is the distance between the expansion points k and l, a,bandcarethecomponentsofthedistance. Charge-chargeinteractionenergy: E ch−ch kl = q k q l R kl (D.1) Charge-dipoleinteractionenergy: E ch−dip kl = q k P x,y,z a μ l a a R 3 kl (D.2) Charge-quadrupoleinteractionenergy: E ch−quad kl = q k P x,y,z a P x,y,z b Θ l ab ab R 5 kl (D.3) Charge-octupoleinteractionenergy: E ch−oct kl = q k P x,y,z a P x,y,z b P x,y,z c Ω l abc abc R 7 kl (D.4) 128 Dipole-dipoleinteractionenergy: E dip−dip kl = P x,y,z a μ k a μ l b R 3 kl −3 P x,y,z a P x,y,z b μ k a μ l b ab R 5 kl (D.5) Dipole-quadrupoleinteractionenergy: E dip−quad kl =−2 P x,y,z a P x,y,z b Θ k ab μ l a b R 5 kl +5 P x,y,z a P x,y,z b Θ k ab ab P x,y,z c μ l c c R 7 kl (D.6) Quadrupole-quadrupoleinteractionenergy: E quad−quad kl =2 P x,y,z a P x,y,z b Θ k ab Θ l ab 3R 5 kl −20 P x,y,z a P x,y,z b Θ k ab b P x,y,z c Θ l ac c 3R 7 kl +35 P x,y,z a P x,y,z b Θ k ab ab 3R 9 kl (D.7) Nuclei-chargeinteractionenergy: E nuc−ch Il = Z I q l r kl (D.8) Nuclei-dipoleinteractionenergy: E nuc−dip Il = Z I P x,y,z a μ l a a R 3 kl (D.9) Nuclei-quadrupoleinteractionenergy: E nuc−quad Il = Z I P x,y,z a P x,y,z b Θ l ab ab R 5 kl (D.10) Nuclei-octupoleinteractionenergy: E nuc−oct Il = Z I P x,y,z a P x,y,z b P x,y,z c Ω l abc abc R 7 kl (D.11) 129 AppendixE Buckinghammultipoles TheconversionfromthespericaltotheBuckinghammul- tipoles The Buckingham charge q in terms of the spherical multipole Q with l = 0 angular momentum: q =Q 0,0 (E.1) The Buckingham dipoleμ in terms of the spherical multipoleQ withl = 1 angular momentum: μ x =−2Q 11 μ y =2Q 1 ¯ 1 μ z =Q 10 (E.2) 130 The Buckingham quadrupole Θ in terms of the spherical multipole Q with l = 2 angularmomentum: Θ xx =6Q 22 −Q 20 Θ xy =−6Q 2 ¯ 2 Θ yy =−6Q 22 −Q 20 Θ xz =−3Q 21 Θ yz =3Q 2 ¯ 1 Θ zz =2Q 20 (E.3) TheBuckinghamoctupoleΩintermsofthesphericalmultipoleQwithl =3angular momentum: Ω xxx =−30Q 33 +6Q 31 Ω xxy =30Q 3 ¯ 3 −2Q 3 ¯ 1 Ω xxz =−3Q 30 +10Q 32 Ω xyy =30Q 33 +6Q 31 Ω xyz =−10Q 3 ¯ 2 Ω xzz =−8Q 31 Ω yyy =−30Q 3 ¯ 3 −6Q 3 ¯ 1 Ω yyz =−3Q 30 −10Q 32 Ω yzz =8Q 3 ¯ 1 Ω zzz =6Q 30 (E.4) 131 The conversion from the Cartesian to the Buckingham multipoles . TheBuckinghamchargeq intermsoftheCartesianchargeq 0 : q =q 0 (E.5) TheBuckinghamdipoleμintermsoftheCartesiandipoleμ 0 : μ x =μ 0 x μ y =μ 0 y μ z =μ 0 z (E.6) TheBuckinghamquadrupoleΘintermsoftheCartesianquadrupoleΘ 0 : halftrace= Θ 0 xx +Θ 0 yy +Θ 0 zz 2 Θ xx = 3 2 Θ 0 xx −halftrace Θ xy = 3 2 Θ 0 xy Θ yy = 3 2 Θ 0 yy −halftrace Θ xz = 3 2 Θ 0 xz Θ yz = 3 2 Θ 0 yz Θ zz = 3 2 Θ 0 zz −halftrace (E.7) 132 TheBuckinghamoctupoleΩintermsoftheCartesianoctupoleΩ 0 : traceX =Ω 0 xxx +Ω 0 xyy +Ω 0 xzz traceY =Ω 0 xxy +Ω 0 yyy +Ω 0 yzz traceZ =Ω 0 xxz +Ω 0 yyz +Ω 0 zzz Ω xxx = 5 2 Ω 0 xxx − 3 2 traceX Ω xxy = 5 2 Ω 0 xxy − 3 2 traceY Ω xxz = 5 2 Ω 0 xxz − 3 2 traceZ Ω xyy = 5 2 Ω 0 xyy − 1 2 traceY Ω xyz = 5 2 Ω 0 xyz − 1 2 traceZ Ω xzz = 5 2 Ω 0 xzz − 1 2 traceX Ω yyy = 5 2 Ω 0 yyy − 1 2 traceX Ω yyz = 5 2 Ω 0 yyz − 1 2 traceZ Ω yzz = 5 2 Ω 0 yzz − 1 2 traceY Ω zzz = 5 2 Ω 0 zzz (E.8) 133 AppendixF TheQ-CheminputforEFPwater cluster The Q-Chem input file a for a molecular system with two effective fragment water moleculesandanotherwatermoleucleinthetheactiveregionispresentedbelow: $molecule 0 1 O 2.98679899864 0.000000000000 0.000000000000 H 3.37280799847 0.408736299796 0.759199899653 H 3.37280799847 0.408736299796 -0.759199899653 $end $rem basis 6-31G * exchange hf $end $efp_fragments water 2.5 0.0 0.0 water 0.0 2.5 0.0 $end $efp_params fragment water o 2.98679899864 0.000000000000 0.000000000000 h 3.37280799847 0.408736299796 0.759199899653 h 3.37280799847 0.408736299796 -0.759199899653 mult 2.98679899864 0.000000000000 0.000000000000 8.00000 mult 3.37280799847 0.408736299796 0.759199899653 1.00000 mult 3.37280799847 0.408736299796 -0.759199899653 a LinesendingatabackslashsymbolshouldbejoinedwiththenextlineintheactualQ-Cheminput. 134 1.00000 mult 2.98679899864 0.000000000000 0.000000000000 cdamp 2.432093890 -5.9880774156 0.9069997308 0.960401737 0.000000000 -2.0776300674 -2.055053110 -1.094067772 0.1972105926\ 0.0000000000 0.000000000 1.2431711201 1.307173512 0.0000000000 0.3882335020\ 0.0000000000 0.357964785 0.0000000000 0.3997408222\ 0.4232766191 0.000000000 mult 3.37280799847 0.408736299796 0.759199899653 cdamp 1.856169014 -0.6161272254 -0.0607832521 -0.0643620268 -0.1385305992 -0.2984353445 -0.2970178207 -0.2629326306 0.0123821246\ 0.0186582649 0.0197568195 -0.1419121249 -0.1496397602 -0.2838114377 -0.0466364802\ -0.1010324765 -0.0434504063 -0.0996406903 -0.0238456713\ -0.0252496482 0.0121573058 mult 3.37280799847 0.408736299796 -0.759199899653 cdamp 1.856168614 -0.6161272254 -0.0607832521 -0.0643620268 0.1385305992 -0.2984353445 -0.2970178207 -0.2629326306 0.0123821246\ -0.0186582649 -0.0197568195 -0.1419121249 -0.1496397602 0.2838114377 -0.0466364802\ 0.1010324765 -0.0434504063 0.0996406903 -0.0238456713\ -0.0252496482 -0.0121573058 mult 3.17980349856 0.204368149898 0.379599949853 cdamp 10.000000000 -1.3898340671 0.0671975530 0.0711539868 -0.0800483392 -0.7200126389 -0.7080787920 -0.3459299814 0.1042426037\ 0.1265126185 0.1339613836 0.0823824078 0.0860368352 -0.1633636547 0.0224998491\ -0.0684206741 0.0201192189 -0.0686914603 0.0517621518\ 0.0548097854 -0.0023653270 mult 3.17980349856 0.204368149898 -0.379599949853 cdamp 10.000000000 -1.3898340671 0.0671975530 0.0711539868 0.0800483392 -0.7200126389 -0.7080787926 -0.3459299814 0.1042426037\ -0.1265126185 -0.1339613836 0.0823824078 0.0860368352 0.1633636547 0.0224998491\ 0.0684206741 0.0201192189 0.0686914603 0.0517621518\ 135 0.0548097854 0.0023653270 pol 3.225043346880 0.252271132154 -0.448809812124 0.3364082841 0.3771885556 2.0857377185 0.3562152772\ -0.6146554330 -0.6508454378 0.3562153645 -0.9006532748\ -0.9536819287 pol 3.22504358327 0.252271266354 0.448809784236 0.3364078945 0.3771881750 2.0857352855 0.3562148807\ 0.6146555083 0.6508456481 0.3562149890 0.9006533941\ 0.9536821303 pol 2.68786044569 0.0680155240902 0.000000000000 0.0546318425 0.2156829505 0.4678560525 0.0987167685\ -0.0000001586 -0.0000002666 0.1193632771 -0.0000002037\ -0.0000002725 pol 3.07178559579 -0.294563387213 0.000000000000 0.2027164107 0.0676015168 0.4678627585 0.1278414717\ 0.0000000833 0.0000000562 0.1071947676 0.0000000844\ 0.0000000709 disp 3.22504334688 0.252271132154 -0.448809812124 0.9330994836 0.9328329368 0.9312726852 0.9258111592\ 0.9104538267 0.8718871966 0.7834519894 0.6086537072\ 0.3529710151 0.1263476070 0.0217816395 0.0007641188 disp 3.22504358327 0.252271266354 0.448809784236 0.9330984158 0.9328318692 0.931271618733 0.92581009\ 0.9104527744 0.8718861726 0.7834510398 0.6086529349\ 0.3529705524 0.1263474433 0.0217816120 0.0007641179 disp 2.68786044569 0.068015524090 0.000000000000 0.2460531659 0.2459405339 0.2452820477 0.2429879390\ 0.2366256365 0.2211913563 0.1883852733 0.1320296167\ 0.0651263370 0.0192960190 0.0029443470 0.0001007762 disp 3.07178559579 -0.294563387213 0.000000000000 0.2460564460 0.2459438126 0.2452853183 0.2429911814\ 0.2366288005 0.2211943284 0.1883878314 0.1320314467\ 0.0651272751 0.0192963144 0.0029443942 0.0001007779 er_basis sto-3g er_wavefunction 1 9.12366435E-02 -2.55851988E-01 -2.74993695E-01 -2.91183082E-01 1 4.26891343E-01 1.09661524E-01 -5.16313349E-01 2 9.12367358E-02 -2.55852205E-01 -2.74994306E-01 -2.91183344E-01 2 -4.26890773E-01 -5.16313220E-01 1.09660818E-01 136 3 1.54680314E-01 -6.48520365E-01 7.87066922E-01 -1.96453849E-01 3 -1.20289533E-07 8.88260577E-02 8.88262341E-02 4 1.54680523E-01 -6.48521682E-01 -2.41107828E-01 7.74554226E-01 4 3.89668406E-08 8.88271814E-02 8.88271243E-02 er_fock_matrix -0.8097112963 -0.1811928224 -0.8097117054 -0.1951745608 -0.1951748560 -0.5598738582 -0.1951742273 -0.1951744336 -0.1696642207 -0.5598736581 er_lmos 3.22504334688 0.252271132154 -0.448809812124 3.22504358327 0.252271266354 0.448809784236 2.68786044569 0.068015524090 0.000000000000 3.07178559579 -0.294563387213 0.000000000000 $end 137 AppendixG Q-Chem’skeywordsaffectingtheEFP parameterscomputation Q-Chem’s keywords affecting the computation of the electrostatics and polarization parametersarepresentedbelowintheformatcompartibliewiththeQ-Chemmanual: POL GEN Whetherdistributedpolarizabilitiesshouldbecomputedformakeefpjobtype. TYPE: LOGICAL DEFAULT: TRUE OPTIONS: TRUE/FALSE RECOMMENDATION: none POL PRINT BOHR PrintpolarizabilitiesinBohrinsteadofAngstroms. TYPE: LOGICAL DEFAULT: FALSE OPTIONS: TRUE/FALSE RECOMMENDATION: none 138 POL DIFF TYPE Finitedifferentitationtype TYPE: INTEGER DEFAULT: 2 OPTIONS: 2 twostepsfinitedifference(forhigheraccuracy) 1 onestepfinitedifference(faster) RECOMMENDATION: none POL FIELD DELTA Electricfiled(stepsizefornumericaldifferentiationofpolarizabilities) TYPE: INTEGER DEFAULT: 100 OPTIONS: n n∗10 −7 stepsize. RECOMMENDATION: none POL LMO FREEZE HowmanycoreorbitalstofreezeforLMOcomputation TYPE: INTEGER DEFAULT: -1 OPTIONS: -1 freezeallcoreorbitals n freezencoreorbitals RECOMMENDATION: none DMA GEN Whethertoperformdistributedmultipoleanalysis(DMA)forjobtypemakeefp. TYPE: LOGICAL DEFAULT: TRUE OPTIONS: TRUE/FALSE RECOMMENDATION: none 139 DMA MIDPOINTS Add bond midpoint to DMA expansion points on bond centers (in addtion to expansion pointsonnucleai) TYPE: LOGICAL DEFAULT: TRUE OPTIONS: TRUE/FALSE RECOMMENDATION: none DAMP STO1G GUESS WhethertouseSTO1Gguessprocedure TYPE: LOGICAL DEFAULT: TRUE OPTIONS: TRUE/FALSE RECOMMENDATION: none DAMP GEN EXP FIT Whethertocomputeexponentialdamping TYPE: LOGICAL DEFAULT: FALSE OPTIONS: TRUE/FALSE RECOMMENDATION: none DAMP GEN GAUSS FIT Whethertocomputegaussiandamping TYPE: LOGICAL DEFAULT: FALSE OPTIONS: TRUE/FALSE RECOMMENDATION: none 140 DAMP EXP TOL Toleranceforexponentialfitting TYPE: INTEGER DEFAULT: 10 OPTIONS: n n∗10 −7 tolerance RECOMMENDATION: none DAMP GAUSS TOL Toleranceforgaussianfitting TYPE: INTEGER DEFAULT: 1000 OPTIONS: n n∗10 −7 tolerance RECOMMENDATION: none DAMP GRID STEP Gridstepfordampingofelectrostaticpotential TYPE: INTEGER DEFAULT: 50 OPTIONS: n n/100Angstroms RECOMMENDATION: none DAMP RMIN Minimumdistancefromanucleaitoagridpoint TYPE: INTEGER DEFAULT: 67 OPTIONS: n n/100Angstroms RECOMMENDATION: none 141 DAMP RMAX Maximumdistancefromanucleaitoagridpoint TYPE: INTEGER DEFAULT: 300 OPTIONS: n n/100Angstroms RECOMMENDATION: none 142
Abstract (if available)
Abstract
Method developments and applications of the condensed phase and gas phase modeling techniques are presented in two parts of this work.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Electronic structure and spectroscopy of radicals in the gas and condensed phases
PDF
Electronic structure and spectroscopy of excited and open-shell species
PDF
Spectroscopic signatures and dynamic consequences of multiple interacting states in molecular systems
PDF
Electronic structure of ionized non-covalent dimers: methods development and applications
PDF
Electronically excited and ionized states in condensed phase: theory and applications
PDF
Computational spectroscopy in gas and condensed phases
PDF
Ultrafast spectroscopy of aromatic amino acids and their chromophores in the condensed phase
PDF
Spectroscopy of hydrogen and water clusters in helium droplets
PDF
Development of predictive electronic structure methods and their application to atmospheric chemistry, combustion, and biologically relevant systems
PDF
Properties of hard-core bosons in potential traps
PDF
Theoretical modeling of nanoscale systems: applications and method development
PDF
Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics
PDF
Electronic structure of strongly correlated systems: method development and applications to molecular magnets
PDF
The synthesis and characterization of [Ga]ZSM-5 by NMR spectroscopy and density functional theory
PDF
Modeling x-ray spectroscopy in condensed phase
PDF
Two-photon absorption spectroscopy and excited state photochemistry of small molecules
PDF
Ab initio methodologies in studying enzymatic reactions
PDF
Phenomenological AdS/CFT: Applications to condensed matter
PDF
Monetary policy and the term structure of interest rates
PDF
Liquid phase photoelectron spectroscopy as a tool to probe excess electrons in liquid ammonia along with Birch chemistry carbanions: from dilute electrolytic solutions to the metallic regime
Asset Metadata
Creator
Vanovschi, Vitalii
(author)
Core Title
Electronic structure and spectroscopy in the gas and condensed phase: methodology and applications
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Chemistry
Publication Date
02/20/2011
Defense Date
01/09/2009
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
benzene clusters,coupled-cluster methods,density functional theory,effective fragment potential,OAI-PMH Harvest,para-benzyne radical anion,photoelectron spectroscopy,symmetry breaking
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Krylov, Anna I. (
committee chair
), Haas, Stephan (
committee member
), Wittig, Curt (
committee member
)
Creator Email
vanovsky@gmail.com,vvanovsc@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m1980
Unique identifier
UC1230125
Identifier
etd-Vanovschi-2613 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-147045 (legacy record id),usctheses-m1980 (legacy record id)
Legacy Identifier
etd-Vanovschi-2613.pdf
Dmrecord
147045
Document Type
Dissertation
Rights
Vanovschi, Vitalii
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
benzene clusters
coupled-cluster methods
density functional theory
effective fragment potential
para-benzyne radical anion
photoelectron spectroscopy
symmetry breaking