Noncovalent interactions between molecules spectra of weakly bound molecular clusters and

Noncovalent interactions between molecules
and
spectra of weakly bound molecular clusters
Ad van der Avoird
Theoretical Chemistry
Outline:
• Preamble: Born-Oppenheimer approximation
• What are non-covalent interactions?
Quantum mechanical derivation
• How to compute intermolecular force fields ab initio?
• How to test intermolecular force fields?
Van der Waals molecules, spectra
• Illustration: ab initio force field for water, applications
• Molecular collisions, scattering resonances
Useful concepts:
• Molecular force fields (MM calculations)
• Interatomic / intermolecular forces (potential energy surfaces)
• Forces on atoms in solids
• Equilibrium structure, force constants of a molecule
• Chemical reaction paths (from QM calculations)
Exist only in the Born-Oppenheimer approximation
Born-Oppenheimer (adiabatic) approximation
Step 1:
Solve electronic Schrödinger equation
He φ(r; R) = Ee(R) φ(r; R)
for nuclei fixed at positions R.
Involves neglect of nuclear kinetic energy Tn.
Step 2:
Solve nuclear Schrödinger equation
[Tn + Ee(R)] χ(R) = Eχ(R)
with potential energy surface Ee(R)
⇒ vibrations, rotations, (phonons, librations),
chemical reaction dynamics, molecular collisions.
Alternative for step 2
Molecular dynamics (MD):
solve nuclear motions classically on potential surface Ee(R).
Derivation of the Born-Oppenheimer approximation
Exact (non-relativistic) Hamiltonian
H = Tn + Te + V (r, R)
with
Te = −
X ~2
i 2m
∇2
i
and
Tn = −
X ~2
A 2MA
∇2
A
X ZAe2
X
e2
ZAZB e2
−
+
V (r, R) =
i>j |ri − rj |
i,A |ri − RA |
A>B |RA − RB |
X
Electronic Hamiltonian (with the nuclei fixed at positions R)
He = Te + V (r, R)
The total Schrödinger equation reads
H Ψ(r, R) = E Ψ(r, R)
Expand the total wave function
Ψ(r, R) =
X
φk (r; R)χk (R)
k
in solutions φk (r; R) of the electronic Schrödinger equation
He φk (r; R) = Ek (R) φk (r; R)
and substitute it into the total Schrödinger equation.
Multiply by the function φk′ (r; R) from the left and integrate over the
electronic coordinates r. The electronic Hamiltonian He is diagonal
hφk′ (r; R)|He |φk (r; R)i(r) = δk′k Ek (R)
and the electronic wave functions are orthogonal
h φk′ (r; R) | φk (r; R)i(r) = δk′k
This yields a set of coupled eigenvalue equations for the nuclear
wave functions
[ Tn + Ek′ (R) − E ] χk′ (R) =
X
[ Fn]k′k χk (R)
k
Coupling between different electronic states k′, k
[ Fn]k′k (R) = hφk′ (r; R)|Tn| φk (r; R)i(r) − Tn δk′k
occurs through the nuclear kinetic energy operator Tn.
When this coupling is neglected one obtains the Born-Oppenheimer
nuclear Schrödinger equation for electronic state k′
[Tn + Ek′ (R)] χk′ (R) = E χk′ (R)
The (non-adiabatic) coupling terms are
[ Fn]k′k (R) = −
X ~2 A 2MA
2hφk′ | ∇Aφk i(r) · ∇A + hφk′ | ∇2
A φk i(r)
They are small because of the large nuclear masses MA in the denominator.
In the first term one may write
h
i
hφk′ | ∇A, He |φk i(r)
hφk′ | ∇Aφk i(r) =
,
Ek (R) − Ek′ (R)
which shows that the coupling is small only when the electronic
energies Ek (R) and Ek′ (R) are well separated. This is normally the
case, and the Born-Oppenheimer approximation holds.
For certain geometries R the energies Ek (R) and Ek′ (R) may be
equal: two (or more) electronic states are degenerate. The BornOppenheimer approximation breaks down.
Breakdown of the Born-Oppenheimer approximation
Examples:
• Open-shell systems: radicals, molecules in excited states
Degeneracies at symmetric structures
⇒ Jahn-Teller, Renner-Teller distortions
(Conical) intersections of different potential surfaces
important in photochemistry
• Metals
Jahn-Teller effect in Benzene+
Degenerate ground state E1g
Sixfold symmetry distorted
by Jahn-Teller coupling with
normal modes of e2g symmetry
π molecular orbitals
Potential surfaces
Two adiabatic potentials
corresponding to the E1g state
as functions of the two ν6
e2g normal mode coordinates
Red circle shows the vibrational zero-point level
⇒ dynamic Jahn-Teller effect
Conical intersection
between
excited S1 and
ground state S0
S1
potentials
of ethene
C2H4
S0
le
ng
torsion a
ngle
g
an
i
gg
wa
Fast non-radiative transition to ground state through
non-adiabatic coupling prevents UV radiation damage in DNA
Metals
Aluminium
20
15
Energy (eV)
10
EF
5
0
−5
0
0.2
0.4
States / eV unit cell
0.6
0.8
For metals
Energy
Interatomic potential
Distance R
⇒ electron-phonon (non-adiabatic) coupling
How were intermolecular forces discovered?
Han-sur-Lesse, December 2012 – p. 1
Equation of state of a gas
Ideal gas
pV = kT
Non-ideal gas (Van der Waals, 1873)
a p + 2 V − b = kT
V
repulsion ⇒ b (eigenvolume)
attraction ⇒ a (reduced pressure)
Han-sur-Lesse, December 2012 – p. 2
Virial expansion (density ρ = 1/V )
2
3
p = kT ρ + B2 (T )ρ + B3 (T )ρ + . . .
with
1
B2 (T ) = −
2
Z
0
∞
∆E(R)
exp −
kT
− 1 4πR2 dR
∆E(R) is the interaction energy between two
atoms/molecules as function of their distance R
Han-sur-Lesse, December 2012 – p. 3
Intermolecular forces
Han-sur-Lesse, December 2012 – p. 4
Intermolecular forces
1909 / 1912
Reinganum, Debye: dipole-dipole, attractive when
orientations are averaged over thermal motion
Han-sur-Lesse, December 2012 – p. 4
Intermolecular forces
1909 / 1912
Reinganum, Debye: dipole-dipole, attractive when
orientations are averaged over thermal motion
1920 / 1921
Debye, Keesom: dipole (quadrupole) - induced dipole
(attractive)
Han-sur-Lesse, December 2012 – p. 4
Intermolecular forces
1909 / 1912
Reinganum, Debye: dipole-dipole, attractive when
orientations are averaged over thermal motion
1920 / 1921
Debye, Keesom: dipole (quadrupole) - induced dipole
(attractive)
1927
Heitler & London: Quantum mechanics (QM)
⇒ covalent bonding for singlet H2 (S = 0)
⇒ exchange repulsion for triplet H2 (S = 1)
Han-sur-Lesse, December 2012 – p. 4
Intermolecular forces
1909 / 1912
Reinganum, Debye: dipole-dipole, attractive when
orientations are averaged over thermal motion
1920 / 1921
Debye, Keesom: dipole (quadrupole) - induced dipole
(attractive)
1927
Heitler & London: Quantum mechanics (QM)
⇒ covalent bonding for singlet H2 (S = 0)
⇒ exchange repulsion for triplet H2 (S = 1)
1927 / 1930
Wang, London: QM ⇒ dispersion forces (attractive)
Han-sur-Lesse, December 2012 – p. 4
QM derivation of intermolecular forces
correspondence with classical electrostatics
Han-sur-Lesse, December 2012 – p. 5
QM derivation of intermolecular forces
correspondence with classical electrostatics
Intermezzo:
(Time-independent) perturbation theory
Han-sur-Lesse, December 2012 – p. 5
Schrödinger equation HΦ = EΦ not exactly solvable.
Perturbation theory ⇒ Approximate solutions Ek and Φk
Han-sur-Lesse, December 2012 – p. 6
Schrödinger equation HΦ = EΦ not exactly solvable.
Perturbation theory ⇒ Approximate solutions Ek and Φk
Find simpler Hamiltonian H (0) for which H (0) Φ(0) = E (0) Φ(0)
(0)
(0)
is solvable, with solutions Ek and Φk
“Perturbation” H (1) = H − H (0)
Write H(λ) = H (0) + λH (1) (switch parameter λ)
λ
0
−→
H (0)
H(λ)
(0)
Ek
(0)
Φk
Ek (λ)
−→
−→
Φk (λ)
−→
1
H
Ek
Φk
Han-sur-Lesse, December 2012 – p. 6
Expand
Ek (λ) =
Φk (λ) =
(0)
(1)
2 (2)
Ek + λEk + λ Ek + . . .
(1)
(0)
2 (2)
Φk + λΦk + λ Φk + . . .
Han-sur-Lesse, December 2012 – p. 7
Expand
Ek (λ) =
Φk (λ) =
(0)
(1)
2 (2)
Ek + λEk + λ Ek + . . .
(1)
(0)
2 (2)
Φk + λΦk + λ Φk + . . .
Substitution into H(λ)Φk (λ) = Ek (λ)Φk (λ) and equating
each power of λ yields, after some manipulations
(1)
Ek
(2)
Ek
=
=h
(0)
Φk
|H
(1)
|
(0)
Φk
i
X h Φ(0) | H (1) | Φ(0) ih Φ(0) | H (1) | Φ(0) i
i
i
k
k
i6=k
(0)
Ek
(0)
− Ei
(0)
Used to calculate perturbation corrections of Ek
Han-sur-Lesse, December 2012 – p. 7
First perturbation correction of
(1)
Φk =
(0)
Φk
X h Φ(0) | H (1) | Φ(0) i
i
k
(0)
i6=k
(0)
Ek − Ei
(0)
Φi
Han-sur-Lesse, December 2012 – p. 8
First perturbation correction of
(1)
Φk =
(0)
Φk
X h Φ(0) | H (1) | Φ(0) i
i
k
(0)
i6=k
(0)
Ek − Ei
(0)
Φi
The second order energy may also be written as
(2)
(0)
(1)
Ek = h Φk | H (1) | Φk i
Han-sur-Lesse, December 2012 – p. 8
Molecule in electric field
External potential V (r) = V (x, y, z)
Han-sur-Lesse, December 2012 – p. 9
Molecule in electric field
External potential V (r) = V (x, y, z)
Particles i with charge qi (nuclei qi = Zi e, electrons qi = −e)
Han-sur-Lesse, December 2012 – p. 9
Molecule in electric field
External potential V (r) = V (x, y, z)
Particles i with charge qi (nuclei qi = Zi e, electrons qi = −e)
Hamiltonian H = H (0) + H (1)
with free molecule Hamiltonian H (0)
and perturbation
H (1) =
n
X
i=1
qi V (ri ) =
n
X
qi V (xi , yi , zi )
i=1
Han-sur-Lesse, December 2012 – p. 9
Multipole (Taylor) expansion
∂V
∂V
∂V
V (x, y, z) = V0 + x
+y
+z
+ ...
∂x 0
∂y 0
∂z 0
with electric field F = (Fx , Fy , Fz ) = −grad V
V (r) = V (x, y, z) = V0 − r · F0 + . . .
Han-sur-Lesse, December 2012 – p. 10
Multipole (Taylor) expansion
∂V
∂V
∂V
V (x, y, z) = V0 + x
+y
+z
+ ...
∂x 0
∂y 0
∂z 0
with electric field F = (Fx , Fy , Fz ) = −grad V
V (r) = V (x, y, z) = V0 − r · F0 + . . .
Perturbation operator
H (1) = qV0 − µ · F0 + . . .
with total charge q =
n
X
i=1
qi and dipole operator µ =
n
X
qi ri
i=1
Han-sur-Lesse, December 2012 – p. 10
First order perturbation energy (for ground state k = 0)
(1)
E0
(0)
(0)
= h Φ0 | H (1) | Φ0 i
(0)
(0)
= h Φ0 | − µ · F0 + . . . | Φ0 i
(0)
(0)
= −h Φ0 | µ | Φ0 i · F0 + . . .
= −h µ i · F0 + . . .
Han-sur-Lesse, December 2012 – p. 11
First order perturbation energy (for ground state k = 0)
(1)
E0
(0)
(0)
= h Φ0 | H (1) | Φ0 i
(0)
(0)
= h Φ0 | − µ · F0 + . . . | Φ0 i
(0)
(0)
= −h Φ0 | µ | Φ0 i · F0 + . . .
= −h µ i · F0 + . . .
Energy of permanent dipole h µ i in field F0 .
Same as classical electrostatics, with dipole h µ i.
Han-sur-Lesse, December 2012 – p. 11
Second order perturbation energy
for neutral molecule (q = 0) and field in z -direction
i.e., F0 = (0, 0, F0 ) and H (1) = −µz F0
(2)
E0
=
X h Φ(0) | H (1) | Φ(0) ih Φ(0) | H (1) | Φ(0) i
0
0
i
i
i6=0
(0)
E0
(0)
− Ei


X h Φ(0) | µz | Φ(0) ih Φ(0) | µz | Φ(0) i
0
0
i
i
 F02
= 
(0)
(0)
E
−
E
i6=0
0
i
Han-sur-Lesse, December 2012 – p. 12
Second order perturbation energy
for neutral molecule (q = 0) and field in z -direction
i.e., F0 = (0, 0, F0 ) and H (1) = −µz F0
(2)
E0
=
X h Φ(0) | H (1) | Φ(0) ih Φ(0) | H (1) | Φ(0) i
0
0
i
i
i6=0
(0)
E0
(0)
− Ei


X h Φ(0) | µz | Φ(0) ih Φ(0) | µz | Φ(0) i
0
0
i
i
 F02
= 
(0)
(0)
E
−
E
i6=0
0
i
Same as classical electrostatics: Epol = −12 αF02 ,
with polarizability
Han-sur-Lesse, December 2012 – p. 12
αzz = 2
X h Φ(0) | µz | Φ(0) ih Φ(0) | µz | Φ(0) i
0
0
i
i
(0)
i6=0
Ei
(0)
− E0
Han-sur-Lesse, December 2012 – p. 13
αzz = 2
X h Φ(0) | µz | Φ(0) ih Φ(0) | µz | Φ(0) i
0
0
i
i
(0)
Ei
i6=0
(0)
− E0
The polarizability αzz can also be obtained from the induced
dipole moment. The total dipole moment is
(0)
(1)
(0)
(1)
h Φ0 + Φ0 | µz | Φ0 + Φ0 i =
h
(0)
Φ0
| µz |
(0)
Φ0
i + 2h
(0)
Φ0
| µz |
(1)
Φ0
i+h
(1)
Φ0
| µz |
(1)
Φ0
i
Han-sur-Lesse, December 2012 – p. 13
αzz = 2
X h Φ(0) | µz | Φ(0) ih Φ(0) | µz | Φ(0) i
0
0
i
i
(0)
Ei
i6=0
(0)
− E0
The polarizability αzz can also be obtained from the induced
dipole moment. The total dipole moment is
(0)
(1)
(0)
(1)
h Φ0 + Φ0 | µz | Φ0 + Φ0 i =
h
(0)
Φ0
| µz |
(0)
Φ0
i + 2h
(0)
Φ0
| µz |
(1)
Φ0
i+h
(1)
Φ0
| µz |
(1)
Φ0
i
The (first order) induced dipole moment µind is the second
term. With the first order wave function
(1)
Φ0 =
X h Φ(0) | H (1) | Φ(0) i
0
i
i6=0
(0)
E0
(0)
− Ei
(0)
Φi
Han-sur-Lesse, December 2012 – p. 13
and H (1) = −µz F0 this yields
(0)
(1)
µind = 2h Φ0 | µz | Φ0 i


X h Φ(0) | µz | Φ(0) ih Φ(0) | µz | Φ(0) i
0
0
i
i
 F0
= 2
(0)
(0)
E
−
E
i6=0
0
i
Han-sur-Lesse, December 2012 – p. 14
and H (1) = −µz F0 this yields
(0)
(1)
µind = 2h Φ0 | µz | Φ0 i


X h Φ(0) | µz | Φ(0) ih Φ(0) | µz | Φ(0) i
0
0
i
i
 F0
= 2
(0)
(0)
E
−
E
i6=0
0
i
As in classical electrostatics: µind = αF0 , with the same
formula for the polarizability αzz as above.
Han-sur-Lesse, December 2012 – p. 14
and H (1) = −µz F0 this yields
(0)
(1)
µind = 2h Φ0 | µz | Φ0 i


X h Φ(0) | µz | Φ(0) ih Φ(0) | µz | Φ(0) i
0
0
i
i
 F0
= 2
(0)
(0)
E
−
E
i6=0
0
i
As in classical electrostatics: µind = αF0 , with the same
formula for the polarizability αzz as above.
For arbitrary molecules the direction of the induced dipole
moment µind is not parallel to F0 . The polarizability α is a
second rank tensor with non-zero elements αxy , etc.
Han-sur-Lesse, December 2012 – p. 14
and H (1) = −µz F0 this yields
(0)
(1)
µind = 2h Φ0 | µz | Φ0 i


X h Φ(0) | µz | Φ(0) ih Φ(0) | µz | Φ(0) i
0
0
i
i
 F0
= 2
(0)
(0)
E
−
E
i6=0
0
i
As in classical electrostatics: µind = αF0 , with the same
formula for the polarizability αzz as above.
For arbitrary molecules the direction of the induced dipole
moment µind is not parallel to F0 . The polarizability α is a
second rank tensor with non-zero elements αxy , etc.
For isotropic systems (atoms, freely rotating molecules)
α is diagonal and αxx = αyy = αzz = α.
Han-sur-Lesse, December 2012 – p. 14
Long range interactions between two molecules
Molecules A and B at distance R with no overlap of their
wave functions. Particles i ∈ A and j ∈ B .
Han-sur-Lesse, December 2012 – p. 15
Long range interactions between two molecules
Molecules A and B at distance R with no overlap of their
wave functions. Particles i ∈ A and j ∈ B .
Hamiltonian H = H (0) + H (1) with free molecule Hamiltonian
H (0) = H A + H B and interaction operator
X X qi qj
H (1) =
rij
i∈A j∈B
Han-sur-Lesse, December 2012 – p. 15
Long range interactions between two molecules
Molecules A and B at distance R with no overlap of their
wave functions. Particles i ∈ A and j ∈ B .
Hamiltonian H = H (0) + H (1) with free molecule Hamiltonian
H (0) = H A + H B and interaction operator
X X qi qj
H (1) =
rij
i∈A j∈B
Same as H (1) =
n
X
qi V (ri ) in previous section with
i=1
molecule A in electric potential
X qj
V (ri ) =
rij
j∈B
of molecule B .
Han-sur-Lesse, December 2012 – p. 15
Multipole expansion of the interaction operator
ri = (xi , yi , zi ),
rij = rj − ri + R
rj = (xj , yj , zj ),
and
R = (0, 0, R)
rij = |rij |
Han-sur-Lesse, December 2012 – p. 16
A double Taylor expansion in (xi , yi , zi ) and (xj , yj , zj ) of
h
i−1/2
1
= (xj − xi )2 + (yj − yi )2 + (zj − zi + R)2
rij
at (xi , yi , zi ) = (0, 0, 0) and (xj , yj , zj ) = (0, 0, 0) yields
zj
xi xj + yi yj − 2zi zj
1
1
zi
= + 2− 2+
+ ...
3
rij
R R
R
R
This expansion converges when |ri | + |rj | < R.
Han-sur-Lesse, December 2012 – p. 17
Substitution into H (1) gives, after some rearrangement
A µB + µA µB − 2µA µB
A qB
AqB
A µB
µ
q
µ
q
x x
y y
z z
z
H (1) =
+ z2 −
+
R
R
R2
R3
X
X
A
B
with the total charges q =
qi q =
qj
i∈A
and the dipole operators µA =
j∈B
X
i∈A
qi ri
X
µB =
qj rj
j∈B
Han-sur-Lesse, December 2012 – p. 18
Substitution into H (1) gives, after some rearrangement
A µB + µA µB − 2µA µB
A qB
AqB
A µB
µ
q
µ
q
x x
y y
z z
z
H (1) =
+ z2 −
+
R
R
R2
R3
X
X
A
B
with the total charges q =
qi q =
qj
i∈A
and the dipole operators µA =
j∈B
X
i∈A
qi ri
X
µB =
qj rj
j∈B
This operator H (1) includes the electrostatic interactions
between the charges and dipole moments of the molecules
A and B . Higher (quadrupole) interactions are neglected.
Han-sur-Lesse, December 2012 – p. 18
Alternative forms of the dipole-dipole interaction operator
are
B + µA µB − 2µA µB
B
A · T · µB
µA
µ
µA · µB − 3µA
µ
µ
x x
y y
z z
z z
=
=
3
3
R
R
R3
with the interaction tensor

 
Txx Txy Txz
1

 
T =  Tyz Tyy Tyz  =  0
Tzx Tzy Tzz
0
0
1
0

0

0 
−2
This tensor can also be expressed in more general
coordinates.
Han-sur-Lesse, December 2012 – p. 19
The solutions of the Schrödinger equations of the free
molecules A and B are
A A
H A ΦA
=
E
k1
k1 Φk1
B B
H B ΦB
=
E
k2
k2 Φk2
and of the unperturbed problem
(0)
(0) (0)
H (0) ΦK = EK ΦK
(0)
(0)
B and eigenvalues E
A + EB
with ΦK = ΦA
Φ
=
E
k1 k2
k1
k2
K
Han-sur-Lesse, December 2012 – p. 20
Proof
(0) (0)
H ΦK
=
=
=
A
B
B
H + H ΦA
Φ
k1 k2
B
A
B B
A A
H Φk1 Φk2 + Φk1 H Φk2
A B
A
B
Ek1 + Ek2 Φk1 Φk2
Han-sur-Lesse, December 2012 – p. 21
Proof
(0) (0)
H ΦK
=
=
=
A
B
B
H + H ΦA
Φ
k1 k2
B
A
B B
A A
H Φk1 Φk2 + Φk1 H Φk2
A B
A
B
Ek1 + Ek2 Φk1 Φk2
Perturbation operator (repeated)
H (1)
B
A µB
A · T · µB
q A q B µA
q
q
µ
z
=
+ z2 −
+
2
R
R
R
R3
Each term factorizes in A and B operators !
Han-sur-Lesse, December 2012 – p. 21
The first order energy is
(1)
(0)
(0)
B
(1)
A B
E0 = h Φ0 | H (1) | Φ0 i = h ΦA
Φ
|
H
|
Φ
0 0
0 Φ0 i
Han-sur-Lesse, December 2012 – p. 22
The first order energy is
(1)
(0)
(0)
B
(1)
A B
E0 = h Φ0 | H (1) | Φ0 i = h ΦA
Φ
|
H
|
Φ
0 0
0 Φ0 i
With the multipole expansion of H (1) one can separate
integration over the coordinates (xi , yi , zi ) and (xj , yj , zj )
of the particles i ∈ A and j ∈ B and obtain
(1)
E0
B
A h µB i
A i · T · h µB i
q A q B h µA
iq
q
h
µ
z
z
=
+
−
+
2
2
R
R
R
R3
Han-sur-Lesse, December 2012 – p. 22
The first order energy is
(1)
(0)
(0)
B
(1)
A B
E0 = h Φ0 | H (1) | Φ0 i = h ΦA
Φ
|
H
|
Φ
0 0
0 Φ0 i
With the multipole expansion of H (1) one can separate
integration over the coordinates (xi , yi , zi ) and (xj , yj , zj )
of the particles i ∈ A and j ∈ B and obtain
(1)
E0
B
A h µB i
A i · T · h µB i
q A q B h µA
iq
q
h
µ
z
z
=
+
−
+
2
2
R
R
R
R3
the same as in classical electrostatics, with the permanent
A | ΦA i and
multipole moments h µA i = h ΦA
|
µ
0
0
B | ΦB i
h µB i = h ΦB
|
µ
0
0
Han-sur-Lesse, December 2012 – p. 22
The second order energy is
(2)
E0
=
X h Φ(0) | H (1) | Φ(0) ih Φ(0) | H (1) | Φ(0) i
0
0
K
K
(0)
K6=0
(0)
E0 − EK
The index K that labels the excited states of the system is a
composite index K = (k1 , k2 ). The summation over K 6= 0
can be split into three sums, with
k1 =
6 0, k2 = 0
k1 = 0, k2 =
6 0
k1 6= 0, k2 =
6 0
Molecule A excited
Molecule B excited
Both molecules excited
Han-sur-Lesse, December 2012 – p. 23
(2)
The first term of E0 is
B | H (1) | ΦA ΦB ih ΦA ΦB | H (1) | ΦA ΦB i
X h ΦA
Φ
0 0
0 0
k1 0
k1 0
k1 6=0
E0A − EkA1
Han-sur-Lesse, December 2012 – p. 24
(2)
The first term of E0 is
B | H (1) | ΦA ΦB ih ΦA ΦB | H (1) | ΦA ΦB i
X h ΦA
Φ
0 0
0 0
k1 0
k1 0
k1 6=0
E0A − EkA1
The operator H(1) is term-by-term factorizable and the
integrals in this expression can be separated. For example
A | µA | ΦA ih ΦB | µB | ΦB i
A µB
h
Φ
µ
z
z
0
0
0
k1
z z
B
A B
h ΦA
Φ
|
|
Φ
Φ
i
=
0 0
k1 0
R3
R3
A | ΦA ih µB i
h ΦA
|
µ
z
z
0
k1
=
R3
Furthermore, one may use the orthogonality relation
A i = 0.
h ΦA
|
Φ
0
k1
Han-sur-Lesse, December 2012 – p. 24
A | ΦA i, with the
The transition dipole moments h ΦA
|
µ
z
0
k1
summation over k1 6= 0, occur in the formula for the
A.
polarizability αzz
Han-sur-Lesse, December 2012 – p. 25
A | ΦA i, with the
The transition dipole moments h ΦA
|
µ
z
0
k1
summation over k1 6= 0, occur in the formula for the
A.
polarizability αzz
If one assumes that the polarizability is isotropic,
A = αA = αA = αA , one finds for the first term
αxx
yy
zz
(2)
E0 (pol. A)
αA (q B )2 2αA q B h µB
z i
= −
+
4
2R
R5
2 + h µB i2 + 4h µB i2 )
αA (h µB
i
x
y
z
−
2R6
Han-sur-Lesse, December 2012 – p. 25
Also this results agree with classical electrostatics. The
electric field of the point charge q B at the center of molecule
A is
B
q
F = (Fx , Fy , Fz ) = 0, 0, − 2
R
and the electric field of the permanent dipole moment
h µB i is
!
B i
h
µ
2h µB
h µB
i
y
z i
x
,−
,
F = −
3
3
R
R
R3
(2)
The second order interaction energy E0 (pol. A) is simply
the polarization energy − 12 αA F 2 of molecule A in the
electric field of the charge and dipole of molecule B .
Han-sur-Lesse, December 2012 – p. 26
Analogously, we find for the second term, which includes
a summation over the excited states k2 of molecule B
(2)
E0 (pol. B)
B
(q A )2 αB 2q A h µA
iα
z
= −
−
2R4
R5
2
A 2
A 2 B
(h µA
x i + h µy i + 4h µz i )α
−
2R6
This is the classical energy of polarization of molecule B
in the field of A.
Han-sur-Lesse, December 2012 – p. 27
The third term contains the summation over the excited
states of both molecules. All interaction terms with the
charges q A and q B cancel, because of the orthogonality
A i = 0. Only the dipole-dipole term of H (1)
relation h ΦA
|
Φ
0
k1
is left and we obtain
(2)
E0 (disp)
=
B | H (1) | ΦA ΦB ih ΦA ΦB | H (1) | ΦA ΦB i
X X h ΦA
Φ
0 0
0 0
k1 k2
k1 k2
k1 6=0 k2 6=0
= −R−6
(E0A − EkA1 ) + (E0B − EkB2 )
A
2
A
A
B
B
B
X X h Φ0 | µ | Φk i · T · h Φ0 | µ | Φk i
1
2
(EkA1 − E0A ) + (EkB2 − E0B )
k1 6=0 k2 6=0
Han-sur-Lesse, December 2012 – p. 28
This term, the dispersion energy, has no classical
equivalent; it is purely quantum mechanical.
It is proportional to R−6 .
Han-sur-Lesse, December 2012 – p. 29
This term, the dispersion energy, has no classical
equivalent; it is purely quantum mechanical.
It is proportional to R−6 .
It can be easily proved that each of the three second order
terms is negative. Therefore, the induction and dispersion
energies are always attractive.
Han-sur-Lesse, December 2012 – p. 29
This term, the dispersion energy, has no classical
equivalent; it is purely quantum mechanical.
It is proportional to R−6 .
It can be easily proved that each of the three second order
terms is negative. Therefore, the induction and dispersion
energies are always attractive.
For neutral, non-polar molecules the charges q A , q B and
permanent dipole moments h µA i, h µB i are zero, and the
dispersion energy is the only second order interaction.
Han-sur-Lesse, December 2012 – p. 29
Terms with higher powers of R−1 occur as well. They
originate from the quadrupole and higher multipole
moments that we neglected.
Han-sur-Lesse, December 2012 – p. 30
Terms with higher powers of R−1 occur as well. They
originate from the quadrupole and higher multipole
moments that we neglected.
An approximate formula, due to London, that is often used
to estimate the dispersion energy is
(2)
E0 (disp)
3αA αB I A I B
≈−
2R6 I A + I B
This formula is found if one assumes that all the excitation
energies EkA1 − E0A and EkB2 − E0B are the same, and are
equal to the ionization energies I A and I B .
Han-sur-Lesse, December 2012 – p. 30
Summary of long range interactions
The interactions between two molecules A and B can be
derived by means of QM perturbation theory.
Han-sur-Lesse, December 2012 – p. 31
Summary of long range interactions
The interactions between two molecules A and B can be
derived by means of QM perturbation theory.
The first order energy equals the classical electrostatic
(Coulomb) interaction energy between the charges and
dipole moments of the molecules. It may be attractive or
repulsive, depending on the (positive or negative) charges
and on the orientations of the dipole moments. The dipolar
terms average out when the dipoles are freely rotating.
Han-sur-Lesse, December 2012 – p. 31
The second order energy consists of three contributions.
The first two terms correspond to the classical polarization
energies of the molecules in each other’s electric fields.
The third term is purely QM. All the three contributions are
attractive. They start with R−4 terms when the molecules
have charges and with R−6 terms when they are neutral.
The dispersion energy, with the leading term proportional
to R−6 , occurs also for neutral molecules with no
permanent dipole moments.
Han-sur-Lesse, December 2012 – p. 32
The second order energy consists of three contributions.
The first two terms correspond to the classical polarization
energies of the molecules in each other’s electric fields.
The third term is purely QM. All the three contributions are
attractive. They start with R−4 terms when the molecules
have charges and with R−6 terms when they are neutral.
The dispersion energy, with the leading term proportional
to R−6 , occurs also for neutral molecules with no
permanent dipole moments.
All of these terms can be calculated when the wave
B and energies E A , E B of the free
functions ΦA
,
Φ
k1
k2
k1
k2
molecules A and B are known, but one should somehow
approximate the infinite summations over excited states
k1 and k2 that occur in the second order expressions.
Han-sur-Lesse, December 2012 – p. 32
Interactions in the overlap region
Heitler and London (Valence Bond) wave functions for H2
1sA (r1 )1sB (r2 ) ± 1sB (r1 )1sA (r2 )
with the plus sign for the singlet spin (S = 0) function
α(1)β(2) − β(1)α(2)
and the minus sign for the triplet spin (S = 1) functions
α(1)α(2)
α(1)β(2) + β(1)α(2)
β(1)β(2)
The total electronic wave function is antisymmetric (Pauli)
Han-sur-Lesse, December 2012 – p. 33
Interaction energy ∆E(R) = EH2 − 2EH
Triplet (S = 1)
Q(R) = “Coulomb integral”
Q
−J
_____
∆E
1 − S2
Q
+J
_____
1 + S2
Singlet (S = 0 )
J(R) = “exchange integral”
S(R) = h 1sA | 1sB i
= overlap integral
R
Han-sur-Lesse, December 2012 – p. 34
Interaction is dominated by the exchange integral J(R),
which is negative, so that the exchange interaction is
attractive (covalent bonding) in the singlet state and
repulsive in the triplet state.
Han-sur-Lesse, December 2012 – p. 35
Interaction is dominated by the exchange integral J(R),
which is negative, so that the exchange interaction is
attractive (covalent bonding) in the singlet state and
repulsive in the triplet state.
For He2 there is only one (singlet) state and the interaction
energy ∆E(R) is purely repulsive: exchange (or Pauli)
repulsion or steric hindrance.
Han-sur-Lesse, December 2012 – p. 35
Molecular orbital picture
Covalent
bonding
Exchange
repulsion
Han-sur-Lesse, December 2012 – p. 36
Molecular orbital picture
Covalent
bonding
Exchange
repulsion
H–H interaction
Exchange
repulsion
He–He interaction
Han-sur-Lesse, December 2012 – p. 36
Most stable molecules are closed-shell systems and the
exchange energy between them is always repulsive. It
depends on the overlap between the wave functions of A
and B and decays exponentially with the distance R.
Han-sur-Lesse, December 2012 – p. 37
Most stable molecules are closed-shell systems and the
exchange energy between them is always repulsive. It
depends on the overlap between the wave functions of A
and B and decays exponentially with the distance R.
In combination with attractive long range interactions
(proportional to R−n ) this gives rise to a minimum in ∆E(R).
This, so-called, non-covalent bonding is much weaker than
covalent bonding, except when A and B are (atomic or
molecular) ions with opposite charges (cf. Na+ Cl− ).
Han-sur-Lesse, December 2012 – p. 37
Most stable molecules are closed-shell systems and the
exchange energy between them is always repulsive. It
depends on the overlap between the wave functions of A
and B and decays exponentially with the distance R.
In combination with attractive long range interactions
(proportional to R−n ) this gives rise to a minimum in ∆E(R).
This, so-called, non-covalent bonding is much weaker than
covalent bonding, except when A and B are (atomic or
molecular) ions with opposite charges (cf. Na+ Cl− ).
Binding (merely by the attractive dispersion energy) is
weakest when both molecules are neutral and non-polar:
pure Van der Waals interactions.
Han-sur-Lesse, December 2012 – p. 37
A special type of interactions between polar molecules is
hydrogen bonding:
X–H· · · Y
exchange
repulsion
Binding energy ∆ E
Van der Waals
dispersion
hydrogen bonding
electrostatic (dipole−dipole, etc.)
induction (dipole−induced dipole, etc.)
dispersion
Distance R
No special (HOMO-LUMO, charge-transfer, or weak
covalent bonding) interactions are needed !
Han-sur-Lesse, December 2012 – p. 38
Exercise:
Compute the equilibrium angles of HF–HF at R = 2.75 Å
and H2 O–H2 O at R = 2.95 Å from the dipolar and
quadrupolar interactions only.
Han-sur-Lesse, December 2012 – p. 39
Non-covalent interactions and hydrogen bonding, in
particular, are very important in biology. Alpha helices and
beta sheets in proteins are stabilized by intra- and
inter-molecular hydrogen bonds, and the double stranded
structure of DNA is held together by hydrogen bonds
between the base pairs.
Han-sur-Lesse, December 2012 – p. 40
Non-covalent interactions and hydrogen bonding, in
particular, are very important in biology. Alpha helices and
beta sheets in proteins are stabilized by intra- and
inter-molecular hydrogen bonds, and the double stranded
structure of DNA is held together by hydrogen bonds
between the base pairs.
It is essential that a hierarchy of interactions exists with
binding energies varying over several orders of magnitude.
Interactions in biological systems must be sufficiently strong
to maintain stable structures, but not so strong that they
prevent rearrangement processes (DNA replication, for
instance).
Han-sur-Lesse, December 2012 – p. 40
Non-covalent force fields computed ab initio
•
Supermolecule calculations
•
Symmetry-adapted perturbation theory (SAPT)
Supermolecule calculations
∆E = EAB − EA − EB
Requirements:
1. Include electron correlation, intra- and inter-molecular
(dispersion energy = intermolecular correlation)
2. Choose good basis, with diffuse orbitals (and “bond functions”)
especially to converge the dispersion energy
3. Size consistency. Currently best method: CCSD(T)
4. Correct for basis set superposition error (BSSE) by computing
EA and EB in dimer basis
Symmetry-adapted perturbation theory (SAPT)
Combine perturbation theory with antisymmetrization A (Pauli) to
include short-range exchange effects.
Advantages:
1. ∆E calculated directly.
2. Contributions (electrostatic, induction, dispersion, exchange)
computed individually. Useful in analytic fits of potential surface.
Advantage of supermolecule method:
Easy, use any black-box molecular electronic structure program
Problems in SAPT:
1. Pauli: AH = HA.
Antisymmetrizer commutes with total Hamiltonian H = H (0) +
H (1), but not with H (0) and H (1) separately.
Has led to different definitions of second (and higher) order
energies.
B not exactly known.
2. Free monomer wavefunctions ΦA
and
Φ
k1
k2
Use Hartree-Fock wave functions and apply double perturbation
theory to include intra-molecular correlation, or use CCSD wave
functions of monomers ⇒ Many-body SAPT.
Program packages:
- SAPT2 for pair potentials
- SAPT3 for 3-body interactions
Most difficult: dispersion interactions
First ab initio calculation of He–He binding curve:
Phys. Rev. Letters, 25 (1970)
– H.F. Schaefer, D.R. McLaughlin, F.E. Harris, and B.J. Alder
page 988: De = 12.0 K
– P. Bertoncini and A. C. Wahl
page 991: De = 11.4 K
Present value:
De = 11.01 K = 7.65 cm−1
≈ 0.1 kJ/mol ≈ 10−3 eV ≈ 3.5 × 10−5 Hartree
Can one use DFT methods?
Reviews by Johnson & DiLabio, Zhao & Truhlar
Many different functionals tested with different basis sets
type
dimer
mean error in De
Van der Waals
Rg2, (CH4)2,
(C2H2)2, (benzene)2
40 - 200 %
dipole-induced dipole
CH4–HF,
H2O–benzene, etc.
15 - 100 %
dipole-dipole
(H2CO)2, etc.
10 - 40 %
hydrogen bonded
(NH3)2, (H2O)2,
(HCOOH)2, etc.
3 - 20%
• Some VdW dimers, (benzene)2 for example, not bound
• B971 best, B3LYP worst
• Often best results without BSSE correction, or smaller basis sets
⇒ Right for wrong reason
Basic problems with DFT
1. Exchange repulsion
– Incorrect asymptotic behavior of one-electron potential:
v(r) → exp(−αr) instead of −1/r
– In intermolecular Coulomb energy no self-term present
self-exchange ⇒ spurious attraction
2. Dispersion
– Intrinsically non-local:
cannot be described by local LDA or semi-local GGA methods
DFT with dispersion explicitly included
vdW-DF: M. Dion, H. Rydberg, E. Schroder, D.C. Langreth, B.I. Lundqvist,
Phys. Rev. Lett. 92 (2004) 246401
Ar–Ar interaction
0
−20
Eint (cm−1)
−40
−60
−80
−100
−120
SAPT(DFT)/PBE0
SAPT(DFT)/B97−2
CCSD(T)/CBS
Benchmark
vdW−DF
−140
−160
3.5
4.0
4.5
5.0
R (Å)
5.5
6.0
6.5
From: R. Podeszwa and K. Szalewicz, Chem. Phys. Lett. 412 (2005) 488
DFT+D methods
• Compute interaction energy with “standard” DFT
• Add atom-atom Cnr−n long range dispersion terms with n =
6, . . .
• Switch between short range correlation contained in DFT and
the long range terms with a smooth (parameterized) switch
function
• Parameters from atomic electron densities, polarizabilities, etc.
(Becke & Johnson, Tkatchenko & Scheffler)
• Empirically optimized parameters (Yang, Grimme)
Alternative: SAPT-DFT
Implemented by G. Jansen et al. (Essen) and K. Szalewicz et al.
(Delaware)
• First order SAPT energy (electrostatic + exchange) computed
with monomer densities and density matrices from Kohn-Sham
DFT
• Second-order SAPT energy (induction, dispersion + exchange)
from (time-dependent) coupled perturbed Kohn-Sham response
functions
Only Hartree-Fock like expressions from Many-Body SAPT needed
⇒ better scaling
Caution !
• SAPT-DFT requires XC potential that is good in inner region
and has correct −1/r behavior for r → ∞
• Coupled time-dependent DFT must be used for (frequencydependent) density-density polarizabilities α(r , r ′, ω)
Both groups, K. Szalewicz (Delaware) and G. Jansen (Essen),
further improved efficiency by implementation of density fitting.
Ar–Ar interaction
From: R. Podeszwa and K. Szalewicz, Chem. Phys. Lett. 412 (2005) 488
0
−20
−1
Eint (cm )
−40
−60
−80
−100
−120
SAPT(DFT)/PBE0
SAPT(DFT)/B97−2
CCSD(T)/CBS
Benchmark
vdW−DF
−140
−160
3.5
4.0
4.5
5.0
R (Å)
5.5
6.0
6.5
A. Heßelmann, G. Jansen, M. Schütz, J. Chem. Phys. 122 (2005) 014103
(benzene)2
1512 GTOs (aug-cc-pVQZ), extrapolation to basis set limit
MP2
−14.4
−15.1
−20.3
CCSD(T)
−6.7
−11.8
−11.4
DF-SAPT-DFT
−7.6
−11.9
−12.7
unbound
metastable
unbound
standard DFT
kJ/mol
Adenine-Thymine (G. Jansen et al.)
DF-SAPT-DFT up to aug-cc-pVQZ level
-100 -50
E [kJ/mol]
0
50
100
150
-100 -50
(1)
150
Eel
(1)
(1)
Eexch
Eexch
(2)
(2)
Eind
Eind
Eexch-ind
Eexch-ind
(2)
(2)
(2)
(2)
Edisp
Edisp
Eexch-disp
Eexch-disp
(2)
Eint
100
(1)
Eel
δ(HF)
E [kJ/mol]
0
50
(2)
DFT-SAPT
MP2
CCSD(T)
δ(HF)
Eint
DFT-SAPT
MP2
CCSD(T)
40
Pair interactions
exch
in water trimer
30
kcal/mole
20
10
disper
0
-10
induct
-20
-30
elec
total
(per hydrogen bond)
Tetramer
kcal/mole
0
0
-1
-1
-1
-2
-2
-2
-3
-3
-3
-4
-4
-4
-5
-5
-5
-6
-6
-6
-7
-7
-7
bo
ir
pa
0
4- dy
bo
5- dy
bo
t o dy
ta
l
Pentamer
3bo
dy
4bo
dy
to
ta
l
pa
ir
l
ta
to
bo
3-
pa
ir
dy
Trimer
3-
Many-body interactions
How to test non-covalent force fields?
Molecular beam spectroscopy of Van der Waals molecules
W.L. Meerts, Molecular and Laser Physics, Nijmegen
Intermolecular potential
⇓
Cluster (Van der Waals molecule)
quantum levels,
i.e.,
eigenstates of
nuclear motion Hamiltonian
⇓
Van der Waals spectra
Nuclear motion Hamiltonian H = T + V for “normal”
(= semi-rigid) molecules
• single equilibrium structure
• small amplitude vibrations
Use rigid rotor/harmonic oscillator model
For (harmonic) vibrations
• Wilson GF -matrix method ⇒ frequencies, normal coordinates
Rigid rotor model ⇒ fine structure (high resolution spectra)
Nuclear motion Hamiltonian H = T + V for weakly bound
complexes (Van der Waals or hydrogen bonded)
• multiple equivalent equilibrium structures
(= global minima in the potential surface V )
• small barriers ⇒ tunneling between minima
• large amplitude (VRT) motions:
vibrations, internal rotations, tunneling
(more or less rigid monomers)
• curvilinear coordinates
⇒ complicated kinetic energy operator T
Method for molecule-molecule dimers
H2O–H2O, NH3–NH3, C6H6–C6H6, etc.
Hamiltonian
h̄2 ∂ 2
(J − jA − jB )2
R+
+ V (R, ωA, ωB )
H = TA + TB −
2µR ∂R2
2µR2
Monomer Hamiltonians (X = A, B):
2
2
TX = AX (jX )2
+
B
j
+
C
j
(
)
(
)
X
X
X
X
a
c
b
Basis for bound level calculations
(J)
χn(R) DM K (α, β, 0)∗
X
(j )
(j )
DmAk (ωA)∗DmB k (ωB )∗h jA, mA; jB , mB | jAB , K i
A A
B B
mA ,mB
Wave function expansion
Ψk (R, ωA, ωB , Θ, Φ) =
X
| i i cik
i
yields matrix eigenvalue problem
Hck = Ek ck
with
Hij = h i | H | j i
dimension ≤ 300 000 for water dimer
Lanczos/Davidson iterative methods
⇒ lowest 20 eigenvalues Ek and eigenvectors ck
(0)
Start with trial vectors xk
(0)
(1)
Calculate Hxk ⇒ xk
·
·
· (n−1)
(n)
Calculate Hxk
⇒ xk
·
·
(n) ·
Iterate until xk converged ⇒ ck , Ek
The permutation-inversion (PI) symmetry group
For semi-rigid molecules
Use Point Group of Equilibrium Geometry to describe the (normal
coordinate) vibrations
N.B. This point group is isomorphic to the P I group, which contains all “feasible” permutations of identical nuclei, combined with
inversion E ∗.
Molecule
Point group
P I group
C2v
{E, E ∗, (12), (12)∗}
C3v
{E, (123), (132),
(12)∗, (13)∗ , (23)∗}
Example: H2O
P I operation
frame rotation
point group operation
(12)
Rz (π) = C2z
C2z
E∗
Ry (π) = C2y
σxz reflection
(12)∗
Rx(π) = C2x
σyz reflection
permutation ⇒ frame rotation + point group rotation
permutation-inversion ⇒ frame rotation + reflection
Hence: P I-group ≃ point group
For “floppy” molecules/complexes
• multiple equivalent minima in V
• low barriers: tunneling between these
minima is “feasible”.
⇒ additional “feasible” P I-operations
Example
NH3
semi-rigid NH3
+ inversion
tunneling
(umbrella
up ↔ down)
Additional feasible P I-operations
P I(C3v ) =
{E, (123), (132),
(12)∗ , (13)∗, (23)∗ }
P I(D3h ) =
P I(C3v ) ⊗ {E, E ∗}
Also (12), (13), (23)
and E ∗ are feasible
⇐⇒
observable tunneling
splittings in spectrum
For H2O–H2O
PI group G16 = {E, P12} ⊗ {E, P34} ⊗ {E, PAB } ⊗ {E, E ∗}
Equilibrium geometry has Cs symmetry
⇒ 8-fold tunneling splitting of rovib levels
Illustration:
Ab initio water potential
– Tested by spectroscopy on dimer and trimer
– Used in MD simulations for liquid water
R. Bukowski, K. Szalewicz, G.C. Groenenboom, and A. van der Avoird,
Science, 315, 1249 (2007)
Polarizable water pair potential: CC-pol
• From CCSD(T) calculations in aug-cc-pVTZ + bond function
basis
• Extrapolated to complete basis set (CBS) limit at MP2 level
• 2510 carefully selected water dimer geometries
• Estimated uncertainty < 0.07 kcal/mol
(same as best single-point calculations published)
CC-pol: Analytic representation
• Site-site model with 8 sites (5 symmetry distinct) per molecule
– Coulomb interaction,
– dispersion interaction,
– exponential ‘overlap’ terms: first-order exchange repulsion,
second-order exchange induction + dispersion
• Extra polarizability site for induction interaction
• Long range R−n contributions computed by perturbation theory,
subtracted before fit of short range terms
• Good compromise between accuracy of reproducing computed
points (rmsd of 0.09 kcal/mol for ∆E < 0) and simplicity needed
for molecular simulations
Water cluster spectra (far-infrared, high-resolution)
from Saykally group (UC Berkeley)
Used for test of potential:
Pair
potential
⇒
Dimer
VRT levels
⇒
Dimer
spectrum
Pair + 3-body
potential
⇒
Trimer
VRT levels
⇒
Trimer
spectrum
Water dimer tunneling pathways
Acceptor Tunneling
Donor-Acceptor Interchange Tunneling
Donor (Bifurcation) Tunneling
−
B2
−
E2
−
A2
Water dimer
tunneling
levels
(J = K = 0)
+
B1
+
E1
+
A1
Acceptor
Interchange
Bifurcation
tunneling
tunneling
tunneling
(H O) levels computed from various potentials (J=K=0)
2
2
20
B−
B−
2
−
2
−
E
A−
E
−
A2
R. S. Fellers et al., J. Chem. Phys. 110, 6306 (1999)
2
15
−
cm−1
B2
−
B2
−
10
E
−
A2
E−
A−2
−
B2
E−
5
A−2
B+1
B+1
+
0
E
A+1
Experiment
ASP−W
E+
B+1
A+1
E
A+1
B+1
+
ASP−S
B+1
E+
A+1
NEMO3
E+
A+1
MCY
Acceptor tunneling
B+,A−
1
+
a(K=0) + a(K=1) = 15.33 (ab initio)
13.92 (experiment)
H2O dimer
tunneling levels
E ,E
0.44
0.41
Rotational constants
A = ∆ − (B+C)/2 = 7.54
7.59
B+C =
E−,E+
0.37
0.37
from CC-pol-8s
potential (2008
22.3
22.3
0.73
0.71
0.66
0.65
and
experiment
B−2
E−
A−2
∆
a(K=1)
0.55
0.54
B+1,A−1
A+1,B−1
B−2,A+2
A−2,B+2
E+,E−
E−,E+
a(K=0)
0.78
0.75
B+1
E+
A+1
J=K=0
A+1,B−1
B−2,A+2
0.414
0.411
Szalewicz et al.)
1
−
J=K=1
J=K=2
A−2,B+2
Acceptor tunneling
a(K=0) a(K=1) a(K=2)
2.04
1.77
D2O dimer
tunneling levels
0.64
0.62
1.48
1.31
(ab initio)
(experiment)
0.028
0.027
Rotational constants
A = ∆ − (B+C)/2 = 4.20
4.17
B+C =
a(K=2)
0.029
0.027
0.361
0.362
from CC-pol-8s
potential
0.039
0.036
and
experiment
∆
0.038
0.036
a(K=0)
0.042
0.039
J=K=0
a(K=1)
B−
2
E−−
A2
B+
1
+
E+
A1
0.035
0.033
J=K=1
B+,A−
1 1
+ −
A1,B1
B−,A+
2 2
− +
A2,B2
E+,E−
−
+
E ,E
J=K=2
+
−
1
1
B1,A1
E++,E−−
A ,B
B−,A+
2 2
E−−,E++
A2,B2
βA 360
270
Acceptor Tunneling
7
+
+
A+
,
B
,
E
1
1
βA 360
270
−3
7
14
14
7
7
14
90
7
0
0
6
3
14
180
7
90
180
270
360
α
−6
βA 360
180
−7
90
180
270
360
α
−
−
A−
,
B
,
E
2
2
−14
14
14
90
7
6
0
0
180
3
7
−6
−3
90
−7
−1
4
270
potential
0
0
90
180
270
360
α
χA
−
A+
/
A
1
2
12
24
105
Donor-Acceptor Interchange Tunneling
12
24
0
12
−5
−5
χA
12
0
105
χA
−5
105
χD
−12
−5
10
105
2
4
−2
−1
20
0
10
−5
0
potential
105
χD
B1+ / B2−
12
0
24
12
0
105
χD
χ
A
E states localized
105
Ex
24
12
−5
−5
χA
0
12
24
0
105
χD
105
χA
24
12
−5
10
24
12
−5
20
0
105
10
−5
0
105
χD
Ey
0
potential
0
105
χD
Donor (Bifurcation) Tunneling
βA 360
6
−3
270
−3
3
6
−3
3
βA 360
12
24
270
180
−6
3
−6
12
180
3
24
90
−3
6
0
0
3
−3
90
6
180
potential
3
90
270
−3
360
βD
0
0
90
180
270
360
βD
+
−
−
+ / E−
A+
/
B
/
A
/
B
/
E
1
1
2
2
Intermolecular
vibrations
Wave functions
270
8
135
4
6
12
−4
−8
−4
βD 180
−6
270
α 360
Acceptor wag (AW)
2
α 360
Acceptor twist (AT)
−1
Donor torsion (DT)
8
90
45
135
180
γA
45
90
β 180
A
135
6
180
3
−3
−3
6
90
0
0
45
90
90
135
α 360
4
8
270
4
−4
90
180
45
90
−4
48
−6
135
180
βA
8
−3
3
−6
−6
45
−8
−3
−6
0
0
DT + AW
18
9
18
α 360
−3
180
γA
Stretch
DT overtone
270
135
−6
90
0
0
−6 9
45
6
−4
8
−8
4
0
0
12
90
−4
90
−8
4
180
4
−4
−8
8
−4
4
4
180
−3
180
γA
0
4
5
6
7
8
R
0
0
−4
4
−8
−4 −8
45
90
135
180
βA
calculated
experiment
calculated
experiment
K=0
K=0
K=1
K=1
160
2
2
H2O dimer
140
intermolecular
AT
120
vibrations
DT2
DT2
AW
1
1
2
1
2
1
AT
AT
2
AW
1
2
2
AW
AW
2
1
2
1
1
100
from CC-pol-8s
cm−1
DT
potential
DT
DT
2
DT
1
1
1
2
1
2
80
and
60
experiment
40
2
2
20
2
GS
0
2
GS
1
GS
1
GS
calculated
experiment
calculated
experiment
K=0
K=0
K=1
K=1
2
DT2
140
2
D2O dimer
120
intermolecular
vibrations
DT2
AT
cm−1
80
1
AT
1
2
2
1
AW
AW
AW
DT
DT
60
1
2
2
1
AT
1
2
2
1
AW
2
1
1
1
potential
1
2
1
100
from CC-pol-8s
DT
DT
2
1
DT
2
1
2
2
and
experiment
40
20
GS
0
GS
2
1
GS
2
1
1
2
GS
1
2
0
−200
−400
Experiment
CC−pol
SAPT−5s
VRT(ASP−W)III
Second virial
coefficient
of water vapor
B(T) [cm3/mol]
−600
−800
−1000
−1200
−1400
−1600
−1800
200
300
400
500
600
T [K]
700
800
900 1000
Water trimer tunneling pathways
Torsional tunneling (flipping)
4
B
3
1
C
6
5
A
2
dud
uud
Bifurcation (donor) tunneling
Water trimer
tunneling
levels
(J = 0)
H2O trimer torsional levels
−1
cm
109.4
108.3
100
k
3
87.1
75.9
75
±2
65.6
56.2
50
47.6
32.5
25
±1
24.6
22.7
19.4
8.4
0
0
Experiment
CC−pol
TTM2.1
VRT(ASP−W)III
D2O trimer torsional levels
120.1
−1
cm
115.1
k
100
106.2
±2
98.1
3
(k=0)
90.3
94.1
81.6
75
64.6
50
3
41.1
35.4
25
0
±2
28.0
±1
8.5
26.3
25.1
9.8
8.0
0
Experiment
CC−pol
TTM2.1
VRT(ASP−W)III
7.3
2.4
MD simulations of liquid water, T = 298 K
3
CC-pol, 2B
CC-pol+NB
CC-pol+NB(ind)
Experiment
2.5
1.6
1.2
g OH
0.8
2
0.4
g OO
r [Angstroms]
0
1
1.5
2
3
4
5
6
7
5
6
7
1.6
1
1.2
g HH
0.8
0.5
0.4
r [Angstroms]
2
3
4
5
r [Angstroms]
0
0
6
7
1
2
3
4
Atom-atom radial distribution functions
Conclusions
• CC-pol pure ab initio
• Predicts dimer spectra better than semi-empirical potentials fitted to these spectra
• Second virial coefficients in excellent agreement with experiment
• CC-pol + 3-body potential gives good trimer spectrum
• Simulations of liquid water with CC-pol + N-body forces predict
the neutron and X-ray diffraction data equally well as the best
empirical potentials fitted to these data
• Important role of many-body forces in liquid water
Nearly tetrahedral coordination:
3.8 hydrogen bonds, only 2.8 with pure pair potential
General conclusion
CC-pol, with the accompanying many-body interaction
model, provides the first water force field that recovers
the dimer, trimer, and liquid properties well
Resonances (quasi-bound states)
in molecular collisions
Why are resonances important?
See: D. W. Chandler, J. Chem. Phys, 132, 110901 (2010)
on shape/orbiting resonances:
– Long life time of collision complex
– Enable recombination reactions A + B → AB
– Probe long range behavior of the potential
Early (1972-1979) observations in scattering for H-Hg by Scoles
et al. and for H-X, H2-X with X = Ar, Kr, Xe by Toennies et al.
(E ≥ 8 cm−1)
Also observed in dimer spectra, e.g. for He-HF by Nesbitt et al.
Anomalous line broadening seen in Ar-CH4 by Roger Miller,
due to resonances, explained by calculations of AvdA et al.
Difficult to observe resonances in molecular (crossed)
beam scattering experiments
• Collision energies too high
• Difficult to scan the collision energy
• Energy spread too large ⇒ too much averaging
New possibilities:
Stark-decelerated molecular beams ⇒ velocity-tuning and state-selection
Gerard Meijer, Berlin
Crossed beam setups
Bas van de Meerakker, Gerard Meijer: Berlin/Nijmegen
Textbook example of scattering resonances
Square well potential
E
0
0
a
E3
E2
E1
−V0
x
Scattering wave function at energy E
ψ =


2 A eiϕ sin(k0 x)
q
for 0 ≤ x ≤ a, with k0 = 2m(E + V0)
√
for a < x,
with k = 2mE

e−ikx + ei(kx+2ϕ)
Match functions and log-derivatives at x = a
Phase shift
ϕ = arctan
"
#
k
π
tan (k0a) −
k0
2
Weight in well region
A2 =
k2
k2 + V0 cos2(k0a)
Life time
1 ∂ϕ
∂ϕ
τ =
=
v ∂k
∂E
Phase shift / π
4
3
2
1
0
Amplitude2 in well
−1
0
1
0.1
0.3
0.2
0.3
0.2
0.3
E / V0
0.5
0
0
0.1
E / V0
10
Life time τ
0.2
5
0
0
0.1
E / V0
Molecule-molecule collisions, theory
• Use same Hamiltonian as in bound state calculation on dimers
• Use same angular basis
• Instead of using radial basis, solve coupled differential equations in
R by propagator method (log-derivative, Numerov, Airy, etc.)
Time independent quantum scattering theory
Channels, molecule-molecule (asymptotic, uncoupled)
| n i = | vAjAmAkA i | vB jB mB kB i
Plane wave in center-of-mass frame
2
~2kn
n ie
, E=
+ ǫn ,
2µ
Expansion in spherical waves for large R
Ψpw
n =|
eiknR =
2π X
ikn LM
L
ik n · R
c)
YLML (R
i(kn R−L π
2)
e
b
kn = knk
−i(kn R−L π
2)
−e
R
b )∗
iLYLML (k
Coupled channels method
c)
Expand wave function in | nLML i = | n iYLML (R
1
ΨnLML =
R
X
n′L′ML′
| n′L′ML′ i Un′L′M ′ ; nLM (R)
L
L
Substitution into Schrödinger equation
−
~2
2µ
R−1
d2
dR2
!
R + ∆H − E ΨnLML = 0
yields matrix problem
U ′′(R) = W (R) U (R)
Coupling matrix
Wn′L′M ′ ; nLML (R) =
L
2µ
′ ′ ′
h
n
L ML |∆H − E| nLML i
2
~
Coupled-channel calculations
• Solve coupled channel equations U ′′(R) = W (R) U (R)
by propagation
• Match U (R) to asymptotic form ⇒ S-matrix
• Calculate cross sections, integral and differential
• S-matrix eigenvalues e2iϕj ⇒ phase shifts ϕj
• Eigenphase sum ϕ =
X
ϕj ,
jumps by π at resonances
j
∂ϕ
1 ∂ϕ
=
• Resonance life time τ =
v ∂k
∂E
Shape / orbiting resonances
Centrifugal barrier
2
V(R) + L(L+1) / (2 µ R )
0
R
V(R)
Feshbach resonances (multichannel)
Atom-diatom (j = diatom rotation, b = rotational constant)
j=1
V(R,θ) + b j(j+1)
0
j=0
R
V(R,θ)
Example 1
OH-X scattering with X = He, Ne, Ar, Kr, Xe
OH-X cross sections
He
relative cross section / %
100
200
300
400
100
80
60
5/2e
40
20
3/2e
5/2f
Ne
relative cross section / %
0
100
80
5/2e
60
40
3/2e
20
5/2f
0
Xe
relative cross section / %
100
80
3/2e
60
40
5/2e
20
5/2f
0
100
200
300
400
mean collision energy / cm-1
More OH-X cross sections
He
relative cross section / %
8
100
200
300
400
2
100
200
300
400
30
relative cross section / %
400
20
1
4
5/2f
10
5/2e
2
1/2f
3/2e
3/2f
0
100
8
2
3/2f
1/2e
6
3/2e
1
4
5/2f
2
5/2e
1/2f
0
0
8
relative cross section / %
300
1/2e
1/2e
6
0
Xe
200
3/2f
0
Ne
100
2
3/2f
6
3/2e
1
4
5/2f
1/2e
5/2e
2
1/2f
0
100
200
300
400
0
100
200
mean collision energy / cm-1
300
400
200
300
400
Example 2
NH3-He scattering
NH3 monomer
Inversion tunneling (umbrella up ↔ down)
semi-rigid NH3
P I(C3v ) symmetry
+ inversion
tunneling
P I(D3h ) symmetry
also E ∗ feasible
Rotation-inversion Hamiltonian
HA = Hrot + Hinv
Hrot = −
~2
2


jx2
Ixx(ρ)
+
jy2
Iyy (ρ)
+
jz2
Izz (ρ)


∂ −1
∂
~2
Hinv = − g(ρ)−1/2
Iρρ (ρ) g(ρ)1/2
+ Vinv (ρ)
2
∂ρ
∂ρ
Moments of inertia
N
ζ = 3m m+
mN
H
2
1
2
2
Ixx(ρ) = Iyy (ρ) = 3mH rNH 2 sin ρ + ζ cos ρ
2
Izz (ρ) = 3mH rNH
sin2 ρ
2
2
2
Iρρ(ρ) = 3mH rNH cos ρ + ζ sin ρ
Metric tensor diagonal, determinant
g(ρ) = Ixx(ρ)Iyy (ρ)Izz (ρ)Iρρ (ρ)
Potential Vinv (ρ) and vibration-inversion levels v2±
−
1
+
1
−
0
+
0
120
NH3
100
levels
(j k± )
Energy (cm−1)
rotationinversion
30−+
30
31−+
31
32−+
32
33−+
33
80
60
20−+
20
40
20
0
10−+
10
00−+
00
ortho
21−+
21
22−+
22
11−+
11
para
Stark effect in para-NH3
Energy
11−
11+
Electric field
| 11− i state (low-field seeking) can be decelerated
⇒ Scattering with velocity-tuned, pure | 11− i para-NH3 or ND3
NH3-He potential surface
• Hodges and Wheatley, J. Chem. Phys. 114, 8836 (2001)
• New potential calculated by CCSD(T)-F12 method in our group
Re (a0)
De (cm−1)
Hodges
6.13
33.45
ours
6.09
35.08
NH3-He potential
Expanded in spherical harmonics
160
140
vLM (R, ρ)YLM (θ, φ)
L,M
Correct R−n dependence for each L
80
−5
100
−10
Θ (Degrees)
vLM (R, ρ) expressed analytically in R, ρ
120
−30
−25
−20
−15
V (R, ρ, θ, φ) =
X
−10
180
60
40
0
−10
20
5
6
7
8
R (a 0 )
9
10
Coupled-channel calculations for NH3-He
• NH3 inversion explicitly included (basis v2 = 0±, 1± )
• Apply full P I(D3h ) symmetry (of collision complex)
para NH3-He has E ′, E ′′ symmetries, parity ±
• Calculate cross sections σjk± ←11− ,
dσ
, phase shifts
dΩ
Elastic and inelastic cross sections, orbiting resonances
3
10
J = 4,6
J=5
Cross sections (A2)
Elastic 11− ← 11−
2
10
1
Inelastic 11+ ← 11−
10
J=4
J=5
J=4
J=5
0
10
−1
10
−4
10
−3
10
−2
10
−1
10
Energy (cm−1)
0
10
1
10
Bound state calculations with Rmax = 20, 25, 30 a0
J=0
J=1
J=2
J=3
J=4
Energy (cm−1)
10
8
6
4
2
0
Diss.
limit
−2
−4
−6
−8
R < 20 R < 25 R < 30
R < 20 R < 25 R < 30
R < 20 R < 25 R < 30
R < 20 R < 25 R < 30
R < 20 R < 25 R < 30
Phase shift, life times
3
2
Phase shift / π
1
0
−1
−2
−3
−4
−5
−6
0.1
J=0
J=1
J=2
J=3
J=4
J=5
J=6
0.2
0.5
Life time τ
20
10
2
5
10
20
5
10
20
−1
40
30
1
Energy (cm )
J=0
J=1
J=2
J=3
J=4
J=5
J=6
0
−10
0.1
0.2
0.5
1
2
−1
Energy (cm )
Our potential
Hodges-Wheatley potential
3
10
3
10
Elastic 11− ← 11−
J=3
Elastic 11− ← 11−
J = 4,6
J=5
2
2
Cross sections (A )
10
1
Inelastic 11+ ← 11−
10
J=4
J=5
J=4
J=5
1
10
Inelastic 11+ ← 11−
0
J = 4, 6
10
0
10
−1
10
J=4
2
Cross sections (A2)
10
−1
−4
10
−3
10
−2
10
−1
10
Energy (cm−1)
0
10
1
10
10
−4
10
−3
10
−2
−1
10
10
−1
Energy (cm )
0
10
1
10
Elastic and inelastic cross sections, Feshbach resonances
400
22− open
22+ open
21− open
21+ open
Elastic 11− ← 11−
Cross sections (A2)
300
200
Inelastic 11+ ← 11− (×150)
100
0
10
20
30
Energy (cm−1)
40
50
Inelastic cross sections, Feshbach and shape resonances
400
22− open
22+ open
21− open
21+ open
300
Cross sections (A2)
Inelastic 21− ← 11− (×80)
Inelastic 22− ← 11− (×80)
200
100
Inelastic 22+ ← 11− (×200)
Inelastic 21+ ← 11− (×200)
0
30
40
Energy (cm−1)
50
Phase shift, life times
5
Phase shift / π
4
22+ channel open
J=0
J=1
J=2
J=3
J=4
J=5
J=6
3
2
1
0
−1
−2
20
21
22
23
24
25
−1
Energy [cm ]
26
27
28
22
23
24
25
−1
Energy [cm ]
26
27
28
Lifetime τ [5.31 ps]
100
J=0
J=1
J=2
J=3
J=4
J=5
J=6
80
60
40
20
0
20
21
Scattering wave functions at E = 24.36 cm−1 (J = 2, E ′ symmetry)
5
Wave function2
2
x 10
11+
22+
22−
1.5
1
0.5
0
5
10
15
20
25
R
200
2 2 + (L=0)
1 1 + (L=1)
2 2 − (L=1)
2 1 + (L=1)
2 2 + (L=2)
1 1 + (L=3)
2 2 − (L=3)
2 2 + (L=4)
100
Wave function
30
0
−100
−200
−300
−400
5
10
15
20
R
25
30
Scattering wave functions at E = 37.28 cm−1 (J = 3, E ′′ symmetry)
5
11+
11−
22+
22−
21+
21−
1.5
Wave function
2
2
x 10
1
0.5
0
5
10
15
100
R (a0)
20
25
2 2 + (L=1)
2 1 − (L=1)
1 1 + (L=2)
2 2 − (L=2)
1 1 − (L=3)
2 1 − (L=3)
1 1 + (L=4)
2 2 − (L=4)
2 1 + (L=4)
2 2 + (L=5)
2 1 − (L=5)
0
Wave function
30
−100
−200
−300
−400
−500
5
10
15
R (a )
0
20
25
30
Conclusions
• Good agreement between theory and velocity-tuned
experiments for OH-X scattering, with X = He, Ne, Ar, Kr, Xe
• Many orbiting and Feshbach resonances found in NH3-He
scattering between 0 and 120 cm−1
• Life times much larger than normal collision times
• Feshbach resonances should be experimentally observable
Acknowledgements
Theory
Krzysztof Szalewicz (Delaware)
Claude Leforestier (Montpellier)
Koos Gubbels, Gerrit Groenenboom (Nijmegen)
Experiments
Rich Saykally (Berkeley)
Ludwig Scharfenberg, Moritz Kirste, Gerard Meijer (Berlin)
Bas van de Meerakker (Nijmegen)
Alexander von Humboldt Stiftung (Research Award 2007)