Jordan

Slide1: 

Potential energy surfaces: the key to structure, dynamics, and thermodynamics K. D. Jordan Department of Chemistry University of Pittsburgh Pittsburgh, PA ACS PRF Summer School on Computation, Simulation, and Theory in Chemistry, Chemical Biology, and Materials Chemistry, June 15-18, 2005

Slide2: 

Jordan Group – May 2005

Slide3: 

Potential energy surfaces (PES) Key to understanding Chemical reactions Dynamics/energy transfer Spectroscopy Thermodynamics Methods of obtaining and representing PES analytical model potentials quantum chemistry (grid of energies) Quantum chemical energies on grid of geometries can be fit to analytical potentials for subsequent use in studies of spectroscopy or dynamics Limited to about 10 atoms “On the fly” methods can handle larger systems

Slide4: 

Example – Lennard-Jones (LJ) clusters R Isomers different minima on potential energy surface number of isomers grows exponentially with # of atoms a and b – permutation-inversion isomers Ea = Eb ≠ Ec Two atoms: Multiple atoms - assume pairwise additive: a b dispersion (van der Waals) repulsion 1 2 3 1 3 2 R E R ε 21/6σ

Slide5: 

Stationary points for all coordinates Xi local minima – curvature positive in all directions 1st order saddle points – curvature – in one direction, + in all others Potential energy surface for a two-dimensional system, i.e., E(x,y) [from Wales] Contour map of PES; M = minimum, TS =1st order saddle point, S = 2nd order saddle point

Slide6: 

Minimization methods Calculus based methods Steepest descent (1st deriv.) only finds “closest” minimum convergence is guaranteed Newton-Raphson (NR) (1st and 2nd deriv.) not guaranteed to converge Quasi-Newton methods (1st and 2nd deriv.) 2nd derivatives can be evaluated numerically by update procedures Eigenmode following (1st and 2nd deriv.) extended range of convergence Monte Carlo (MC) based methods Simulated annealing Start at high T, and gradually lower T Basin-hopping (a hybrid MC/calculus method) Neural network approaches

Slide7: 

Locating the global minimum – major challenge even small clusters can have over 1010 minima! Brute force approaches, e.g., starting from many initial structures, work for only the simplest systems Monte Carlo methods such as basin hopping useful for systems containing 100 or so atoms (very computationally demanding) Easy to find global minimum Hard to find global minimum E E E(kJ/mol) E(kJ/mol) Figures from Energy Landscapes, by D. Wales.

Slide8: 

folded unfolded partially folded Even though my examples are drawn from cluster systems, the issues considered are relevant for a wide range of other chemical and biological systems, e.g., to the “protein folding” problem. The above figure is from Brooks et al., Science (2001). Entropy Protein folding

Slide9: 

Locating transition states and reaction pathways Harder than locating local minima Elastic band and other 1st derivative (gradient)-based methods Eigenmode following (EF) (1st and 2nd deriv). Methods using analytical Hessian (d2E/dxidxj matrix) Methods with approximate Hessian (update methods) EF method

Slide10: 

Disconnectivity diagram Ar13 (from D. Wales) Disconnectivity diagram Ar38 (from D. Wales) Icosahedral FCC Icosahedral Energy (kJ/mol) Energy (kJ/mol)

Slide11: 

Thermodynamics of clusters from Monte Carlo (or MD) simulations Potential energy vs. T, LJ38 C vs. T, (H2O)8 (Tharrington and Jordan) C vs. T, LJ38 (Liu and Jordan) solid liquid FCC Icosahedral C C

Slide12: 

Magic number clusters arrangements of atoms that are especially stable Often connected with high symmetry illustrate several of the issues discussed thus far Mass spectrum of Cn+: magic # at n = 60 (from Kroto) Mass spectrum of (H2O)nH+: magic # at n = 21(from Castleman + Bowen) 21 60 60 6

Slide13: 

Bimodal potential energy distribution Only low-energy cubic species populated at low T Many inherent (non-cubic) structures populated at high T System shuttles back and forth between “solid” (cubic) and “liquid” (non-cubic) structures Pot. Energy distribution for (H2O)8, T ≈ Tmax Densities of local minima of (H2O)n clusters

Slide14: 

IR spectra of (H2O)nH+, n = 2-11, from Duncan, et al., Science, in press Mass spectra alone tell us very little about the structures. Recently, the combination of new experimental techniques plus electronic structure calculations have enabled researchers to establish the structures of many cluster systems. Our own work has focused on H+(H2O)n and (H2O)n- clusters.

Slide15: 

One of the biggest challenges in theoretical/computational chemistry is choosing the suitable approach Model potentials vs. quantum chemistry (each of these has several variants)? Do we need to allow for temperature? Is the dynamics well described classically, or is a quantum treatment required? In modeling vibrational spectra, does the harmonic approximation suffice? Approach to be adopted dictated by the nature of the problem being studied This will be illustrated by considering the protonated water clusters

Slide16: 

Approaches for modeling model potentials (molecular mechanics/force fields) applicable to thousands of atoms generally neglect polarization and not suitable for cases with rearrangement of electrons quantum chemistry tens – few hundred atoms Wavefunction-based vs. DFT QM/MM methods primary region – treated quantum mechanically Secondary region – treated with a force field

Slide17: 

Choice of theoretical approaches for our studies of H+(H2O)n there is no model potential that provides a near quantitative description of the interactions in protonated water clusters → must use quantum chemical methods (DFT or MP2) for the n = 5 - 8 clusters, the dominant species are not the global minima → must include vibrational ZPE and allow for finite T effects → must employ a scheme which can locate all the low-energy minima (not just those we anticipate) for addressing some aspects of the vibrational spectra, it is necessary to go beyond the harmonic approximation

Slide18: 

Quantum Chemistry (electronic structure methods) Hψ = Eψ H = Hamiltonian : contains kinetic energy operator, el.-nuclear interactions, el.-el. Interactions A complicated partial differential equation In general – must introduce approximations Orders of magnitude more expensive than using model potentials Even fastest methods scale as N3, where N = number of atoms Research underway to get O(N) scaling for large systems But not subject to limitations of model potentials Includes polarization Applies to all bonding situations All properties accessible Software: both commercial and public domain programs GAMESS, Spartan, Gaussian 03, NWChem, Jaguar, and many others

Slide19: 

Properties: charge distributions, dipole moments electrostatic potentials polarizabilities geometries – minima and transition states vibrational spectra electronic excitation and photoelectron spectra NMR shifts thermochemistry For complex systems, the other major challenge is the exploration of configuration space Even if one or two isomers dominate under experimental conditions, it may be necessary to examine a very large number of isomers in the electronic structure calculations Accounting for finite T/energy effects

Slide20: 

Structures responsible for observed spectra For the n = 5 - 8 clusters, these are not the global minimum isomers. H+(H2O)2 H+(H2O)4 H+(H2O)5 H+(H2O)6 H+(H2O)8 H+(H2O)3

Slide21: 

Eel (T=0) Eel(T=0)+ ZPE E(T = T’) H(T=T’) G(T=T’) Account for vibrational zero-point energy From electronic structure calculations Population of excited vibrational, rotational levels Account for PΔV = ΔnRT (ideal gas) Include entropy Accounting for finite temperature on cluster stability Optimize geometries Calculate harmonic frequencies

Slide22: 

(H2O)6H+ Isomers with dangling water molecules (low frequencies) favored by ZPE and by entropy Zundel-type ion dominates under the experimental conditions, T  150 K.

Slide23: 

Comparison of calculated and measured vibrational spectra of H+(H2O)6 Excellent agreement between theory and experiment, except that the harmonic, T = 0 K calculations cannot account for the broadening of the OH stretch spectra of H-bonded OH groups. need to account for vibrational anharmonicity (e.g., stretch/bend coupling) probably also need to account for finite T effects on the spectra Expt. Intensity Intensity Theory

Slide24: 

vibrational spectra of H+(H2O)n, n = 6-27 free-OH region of spectra reflect structural transitions at n = 12 and n = 21(Shin et al., Science, 2004) Collapse to a single line in the free OH stretch region

Slide25: 

Lowest-energy n=21 structure found in ab initio geometry optimizations Dodecahedron with H3O+ on surface (blue) and H2O (purple) inside cage 4 H-bonds with interior H2O causes a rearrangment of the H-bonding in the dodecahedron there are only 9 free-OH groups (Castleman's experiments suggested 10) all free-OH associated with AAD waters - explains single lines in free OH stretch If the excess proton placed on interior water, it rapidly jumps to surface.

Slide26: 

Interplay between spectroscopy and dynamics concentration of ions so low cannot obtain spectra by simple absorption Obtain spectra instead by dissociation Calculated vs. expt. spectra of magic # cluster. No transitions observed in H3O+ OH stretch region Predissociation spectroscopy H+(H2O)n H+(H2O)n-1 + H2O

Slide27: 

If the ion does not fall apart on the timescale of the experiment, no signal will be observed. These problems illustrate the interplay between structure, spectra, and dynamics inherent in much of today’s research Cold clusters Spectra dominated by 2-photon absorption Is it possible that H3O+ OH stretch vibrations undergo appreciable shifts with > T? If so, this could turn off the 2-photon absorption. 130 150 170 190 210 T(K) Tm with Ar without Ar free OH Eigen OH 10-6 s. 10-2 s. τ

Slide28: 

Vibrational anharmonicity Diatomic molecule: V(x) = aox2( 1 + a1x3 + a2x4 + …) harmonic anharmonicity E(v) = 1/2 hωe(v+1/2) – ωexe(v+1/2)2 + ωeye(v+1/2)3 + … Polyatomic molecules: diagonal anharmonicity: Viii, Viiii off-diagonal anharmonicity: Viij, Vijk, Viijj. etc. - couple modes x=(R-Re)/Re ωe = harmonic frequency ωexe, ωeye = first two anharmonicity constants Be = rotational constant αe = vibr.-rot. coupling ωe = sqrt(4ao*Be) αe = (a1 + 1)(6Be2/ ωe) ωexe = (5a12/4 – a2)(3Be/2) Dunham expansion: unique mapping between 1D potential and the spectroscopic parameters This mapping is lost for polyatomic molecules Depends on 3rd and 4th derivatives Several transitions of the H+(H2O)n clusters are not well described in the harmonic approximation

Slide29: 

2nd-order vibrational perturbation theory Requires Viij, Vijk, Viiii, Viijj can be calculated with standard electronic structure codes Can’t handle shared proton in H5O2+ x4 term dominates: PT fails Can’t handle “progressions” as in CH3NO2-(H2O) Vibrational SCF (VSCF) can be done using ab initio PES (grids) can’t handle progressions Vibrational CI need a representation of the PES limited to about 12 degrees of freedom Diffusion Monte Carlo methods difficulty in handling excited states Approaches for treating anharmonicity

Slide30: 

CH3NO2-(H2O) – an example of important off-diagonal vibrational anharmonicity Experimental spectrum displays 5 ( 90 cm-1 spacing) transitions in the OH stretch region – only two lines expected This is a consequence of strong OH stretch/water rock coupling Key coupling term: VSAR = kASRQSQAQS Configuration interaction with Hamiltonian including this cubic term and with product basis set A, AR, AR2, S, SR, SR2, etc, accounts for observed spectrum (S = symmetric OH stretch, A= asymm. OH stretch, R = water rock) Note how this coupling results in a band with overall width of several hundred cm-1 Such couplings important for energy redistribution expt. theory-harmonic OH stretch CH stretch From Johnson, Sibert, Jordan and Myshakin, 2004 theory - anharmonic

Slide31: 

(H2O)2 – an example illustrating the importance of vibrational anharmonicity of frequencies, ZPE, geometry acceptor donor donor Intermolecular vibrations Frequencies calculated using the MP2 method. Anharmonicities calculated using 2nd order vibrational PT. Excellent agreement between the calculated anh. frequencies and experiment.

Slide32: 

Changes in bond lengths of (H2O)2 upon vibrationally averaging R E Re Ro Actually, this raises an interesting question concerning the development of model potentials for classical MC or MD simulations. Namely, should one design the potential to give the correct Re or Ro values?

Slide33: 

Various issues concerning electronic structure calculations

Slide34: 

Challenges facing electronic structure theory There is still no reliable method for calculating accurate interaction energies between molecules and extended systems. Example – coronene (7 fused benzene rings) standard QC methods need flexible basis sets to treat dispersion Near linear dependency, large BSSE with basis sets such as aug-cc-pVTZ not clear MP2 is suitable for this problem DFT methods Could use with plane waves (to solve linear dependency and BSSE problems) But inappropriate due to neglect of dispersion DMC would need to run very long to reduce statistical error below a few tenths of a kcal/mol Excess electron in bulk water or even in a (H2O)20 cluster Need very large basis sets and inclusion of high-order correlation effects Solution in this case possible by use of quantum Drude oscillators

Slide35: 

Some considerations concerning model potentials For simulations of large systems, model potentials are essential Typically, these model potentials include Bond-stretch, bend, torsional contributions. Electrostatics (generally using point charges) Pose special challenges for extended or periodic systems Lennard Jones (dispersion plus short-range repulsion) Growing realization that dipole polarizability is important Can greatly increase the cost of the simulations Many of the issues can be illustrated by considerations of models for water.

Slide36: 

Water models TIP3P – 3 atom-centered charges + OO LJ int. TIP4P – 3 charges (-2q displaced from O), + OO LJ int. Dang-Chang (DC) – like TIP4P, but with polarizable center added to M site (0.215 Å from O atom) TTM – 3 charges (-2q at M site), 12-10-6 (AR-12 + BR-10 + CR-6) OO interaction, 3 polarizable sites AMOEBA – atom-centered charges, dipoles, quadrupoles, OO, HH, and OH LJ, 3 polarizable sites Water dimer: interaction energies (kcal/mol) +q +q M, -2q

Slide37: 

In-plane electrostatic potential of the water monomer from MP2 ab initio calculations from and from the DC water model. Distances in Å. Outer contour = 0.005 au = 3 kcal/mol MP2 – in-plane DC model – in plane -0.005 0.005 -0.005 0.005 DC model: q = +0.519 H atoms, -1.038 M site, 0.215 Å from the O atom. M

Slide38: 

In-plane electrostatic potential: DC – MP2. Outer blue contour -0.0005 au = 0.3 kcal/mol. Distances in Å. Perp.-to-plane electrostatic potential: DC – MP2. Outer black contour 0.0005 au = 0.3 kcal/mol. Distances in Å. A three-point charge model cannot realistically describe the electrostatic potential potential of water!! Yet, nearly all simulations of water, ice, and biomolecules in water use models with simple point charge representations of the charge distribution. In these figures the part of the electrostatic potential near the atoms has been cut out.

Slide39: 

Differences between the electrostatic potentials from a distributed multipole analysis with moments through the quadrupole on each atom and from MP2 level calculations. Overall the agreement is excellent except for short distances. In-plane Perp. to plane GDMA-MP2 0 0

Slide40: 

In-plane electrostatic potential: Amoeba – MP2. Outer blue contour -0.0005 au = 0.3 kcal/mol. Distances in Å. Perp.-to-plane electrostatic potential: Amoeba – MP2. Outer light blue contour 0.0005 au = 0.3 kcal/mol. Distances in Å. Amoeba-MP2 Amoeba should give results identical to GDMA. Differences due to change in HOH angle and scaling of the atomic quadrupoles. 0 0

Slide41: 

More on polarization interactions 2-body interactions – interaction between each pair uninfluenced by other molecules Many-body interactions – Interaction between A and B alters interactions between A and C and B and C. A B C Inert gas clusters – many-body effects dominated by dispersion Water clusters – many-body effects dominated by polarization E = E1 + E2 + E3 + … + En In general the series converges rapidly Water clusters – 3-body contributions represent 20 – 30% of the net binding energy Isolated water monomer – dipole moment = 1.85 D Water molecule in liquid water – dipole moment ~ 2.6 D + - + + - - + - + - μAB μBA μij – dipole induced on i by charges on j μAB in turn induces a dipole moment on B. Infinite series!

Slide42: 

Effective 2-body potentials for water, e.g. TIP4P and SPC/E, have charges that give a dipole significantly larger than experiment for the monomer account in an effective mater for polarization effects in bulk water overestimate dipoles of water molecules at interfaces and in clusters Many strategies have been introduced for treating polarization point polarizable sites – induced dipoles fluctuating charges (in-plane polarization only) Drude oscillators – two fictitious charges coupled harmonically If atom-centered polarizable sites are employed, it is essential to damp the short range interactions to avoid unphysical behavior at short distances

Slide43: 

The orbital picture reconsidered. One of the most extensive concepts in chemistry is the orbital picture. This is so deeply engrained that we sometimes forget that for many electron systems orbitals are a construct (result from assuming separability of the wavefunction) In much of chemistry the orbitals that we consider are valence-like These are precisely the orbitals that can be calculated using electronic structure codes and minimal basis sets. H2: bonding σg and antibonding σu Ethylene: bonding π and σ and antibonding π* and σ* In dealing with the spectroscopy of molecules there are also excited states resulting from promoting electrons into Rydberg orbitals These arise from higher energy atomic orbitals and tend to be spatially extended. Rydberg states are very sensitive to the environment of a molecule and may vanish in the condensed phase (recall properties of the particle in the box)

Slide44: 

Excited states HF, H2O, NH3, and CH4 do not display singlet excited states with valence character The valence states “dissolve” in the Rydberg sea (quote from Robin) HCl, H2S, PH3, and SiH4 do display singlet excited states with valence character With the longer XH bonds of the latter, the empty unfilled valence orbitals drop below the Rydberg orbitals and are observed Anions If the anion lies energetically above the neutral (negative electron affinity), the anion lies in the continuum of the neutral plus a free electron This is the case for Be, N2, ethylene, benzene, CH3Cl, etc. Typically the electron falls off (autoionizes) in 10-14 sec. Poses a special challenge for theory Issues connected with unfilled orbitals

Slide45: 

Potential energy curves of CH3Cl and CH3Cl- Decay processes electron detachment dissociation (CH3 + Cl-) 1,1-dichlorethane electron transmission spectrum of – two peaks due to the two σ* orbitals dissociative attachment – one peak due to the lower-lying anion electron attachment from upper anion to fast to give Cl- (results from P. Burrow, Univ. Nebraska)

Slide46: 

Vibrational excitation cross sections for two vibrations of CH3Cl. The peaks are due to resonances (temporary anion states). From P. Burrow. Temporary anions pose a significant challenge to theory Standard variational approaches → collapse onto continuum Several methods have been developed for treating such species The resonance energy is actually complex Eres = Er –i/2Γ Er = resonance position, Γ = width Time dependence exp(-iE*t): complex energy – decays in time

Slide47: 

Electrons bound in electrostatic potentials Most famous case: dipole bound anions An excess electron bound to a (H2O)6 chain The electron is so extended, that it should be possible to develop a one-electron model approach Important interaction terms Exchange/repulsion Polarization (e--water, water-water) Electrostatics [e- - permanent charges on (H2O)] Dispersion – left out of all earlier model potential studies Cannot simply add a C/R6 term, due to extended nature of excess electron. We have developed a Drude model of excess-electron molecule interactions.

Slide48: 

Drude model +q -q charges +q, -q coupled through a force constant k R The position of the -q charge is kept fixed. In the presence of a field, the system has a polarizability of q2/k. An electron couples to the Drude oscillator via qr∙R/r3 ,

Slide49: 

Drude model based on the Dang-Chang water model M site: 0.215 Å from O atom. Negative charge (-1.038e) plus Drude oscillator with q2/k = α = 1.444 Å3 H charge = 0.519e Determined using procedure of Schnitker and Rossky Scaled so that model potential KT energy reproduces ab initio KT result for (H2O)2- b Damping coefficient scaled so that model potential CI energy reproduces ab initio CCSD(T) result for (H2O)2- r - position of electron R - displacement of the Drude oscillator

Single Drude Oscillator:: 

Single Drude Oscillator: Electron orbitals described in terms of s, p Gaussians. { } 3D harmonic oscillator functions { } Wavefunction: in “MO” basis set Multiple Drude Oscillators: Basis set:

Slide51: 

Several strategies for solving fully self-consistent treatment of e--water polarization, e--water dispersion, intramolecular induction, intramolecular dispersion. self-consistent treatment of e--water polarization, e--water dispersion, intramolecular induction. Treat intramolecular dispersion through R-6 terms. self consistent treatment of e--water polarization, e--water dispersion.. Treat intramolecular induction using classical Drude oscillators and treat intramolecular dispersion through R-6 terms.

Slide52: 

Surface state and interior electron bound states of (H2O)20- Considerable interest in these species in light of recent work from the Neumark and Zewail groups. Geometries provided by M. Head_Gordon. The anion is not bound in the KT and Hartree-Fock approximations. Electron binding is a result of correlation effects which cause a large contraction of the excess electron