composed of the sum of two opposing Morse potentials, each with D = 600

kJ/mol and a harmonic vibration frequency of 100 ps’1 (3336 cm’1 ), one

with minimum at 0.1 nm and the other at 0.2 nm. The lowest two levels

are adiabatic tunneling states, which di¬er only slightly in energy (0.630

kJ/mol) with wave functions that are essentially the sum and di¬erence of

the diabatic ground state wave functions of the individual wells. In each of

the adiabatic states the proton actually oscillates between the two wells with

a frequency of (E1 ’ E0 )/h = 1.579 THz; it tunnels through the barrier. The

excited states all lie above the energy barrier; the proton then has enough

energy to move over the barrier.

4.3.2 Expansion on a basis set

Another class of solutions is found by expanding the unknown solution in a

¬nite number of properly chosen basis functions φn :

3 See, e.g., Mavri and Grdadolnik (2001), who ¬t two Morse potentials plus modifying terms to

high-level quantum calculations in the case of acetyl acetone.

4.3 The few-particle problem 83

V,E (kJ/mol)

70

60

50

40

30

20

10

0

0.1 0.125 0.15 0.175 0.2

xproton (nm)

Figure 4.1 Lowest six proton quantum states in a double-well potential (thick

curve), typical for a proton in a symmetric hydrogen bond. Levels and wave

functions were generated with the shooting method (see text). Wave functions

are indicated by alternating solid and dotted thin lines. The energies of the states

are (in kJ/mol): 0:10.504, 1:11.135, 2:25.102, 3:32.659, 4:45.008, 5:58.804.

ψ(x) = cn φn (x). (4.5)

n

The time-independent Schr¨dinger equation now becomes an eigenvalue

o

equation (see Chapter ??):

Hc = »Sc, (4.6)

where S is the overlap matrix

φ— φm dx.

Snm = (4.7)

n

For an orthogonal basis set the overlap matrix is the unit matrix. The

eigenvalue equation can be solved by standard methods (see, e.g., Press et

al., 1992). It is most convenient to diagonalize the basis set ¬rst,4 which

4 Diagonalization is not a unique procedure, as there are more unknown mixing coe¬cients than

equations. For example, mixing two functions requires four coe¬cients, while there are three

conditions: two normalizations and one zero overlap. One unique method is to start with

84 Quantum chemistry: solving the time-independent Schr¨dinger equation

o

is a one-time operation after which the eigenvalue problems become much

easier. Such methods do yield excited states as well, and are extendable

to higher dimensions, but they are never exact: their accuracy depends on

the suitability of the chosen basis set. Still, the solution of many-electron

problems as in quantum chemistry depends on this approach.

An example of the use of a basis set can be found in the study of proton

transfer over hydrogen bonds (Berendsen and Mavri, 1997). For solving

the time-dependent SE (see Chapter 5) a description of the proton wave

function for a ¬‚uctuating double-well potential in terms of a simple basis set

was needed. The simplest basis set that can describe the tunneling process

consists of two Gaussian functions, each resembling the diabatic ground state

solution of one well (Mavri et al., 1993). Many analytical theories rely on

such two-state models. It turns out, however, that a reasonable accuracy can

only be obtained with more than two Gaussians; ¬ve Gaussians, if properly

chosen, can reproduce the ¬rst few eigenstates reasonably well (Mavri and

Berendsen, 1995).

4.3.3 Variational Monte Carlo methods

The variational method consists of de¬ning some ¬‚exible function with a

number of adjustable parameters that is expected to encompass a good

approximation to the true wave function. The variational principle says that

the expectation of the Hamiltonian over any function ψ(r) (which must be

quadratically integrable) is always larger than or equal to the ground state

energy with equality only when the function is identical to the ground state

eigenfunction:

ˆ ˆ

ψ(r)Hψ(r) dr Hψ

≥ E0

E= = (4.8)

ψ 2 (r) ψ

ψ2

ˆ

Therefore the parameter values that minimize the expectation of H yield the

best approximation to the ground state wave function. For low-dimensional

problems and for linear dependence on the parameters, the integral can in

general be evaluated, and the minimization achieved. However, for multidi-

mensional cases and when the integral cannot be split up into a linear com-

bination of computable components, the multidimensional integral is better

one normalized eigenfunction, say φ1 , then make φ2 orthogonal to φ1 by mixing the right

amount of φ1 into it, then normalizing φ2 , and proceeding with φ3 in a similar way, making it

orthogonal to both φ1 and φ2 , etc. This procedure is called Schmidt orthonormalization, see,

e.g., Kyrala (1967). A more symmetrical result is obtained by diagonalizing the overlap matrix

of normalized functions by an orthogonal transformation.

4.3 The few-particle problem 85

solved by Monte Carlo techniques. When an ensemble of con¬gurations is

generated that is representative for the probability distribution ψ 2 , the in-

ˆ

tegral is approximated by the ensemble average of Hψ/ψ, which is a local

property that is usually easy to determine. The parameters are then var-

ied to minimize E . The trial wave function may contain electron“electron

correlation terms, for example in the Jastrow form of pair correlations, and

preferably also three-body correlations, while it must ful¬ll the parity re-

quirements for the particles studied.

The generation of con¬gurations can be done as follows. Assume a starting

con¬guration r 1 with ψ 2 (r 1 ) = P1 is available. The “local energy” for the

starting con¬guration is

ˆ

Hψ(r 1 )

µ1 = . (4.9)

ψ(r 1 )

(i) Displace either one coordinate at the time, or all coordinates simul-

taneously, with a random displacement homogeneously distributed

over a given symmetric range (’a, a). The new con¬guration is r 2 .

(ii) Compute P2 = ψ 2 (r 2 ).

(iii) If P2 ≥ P1 , accept the new con¬guration; if P2 < P1 , accept the new

con¬guration with probability P2 /P1 . This can be done by choosing

a random number · between 0 and 1 and accept when · ¤ P2 /P1 .

(iv) If the move is accepted, compute the “local energy”

¯

Hψ(r 2 )

µ2 = ; (4.10)

ψ(r 2 )

if the move is rejected, count the con¬guration r 1 and its energy µ1

again;

(v) Repeat steps (i)“(iv) N times.

(vi) The expectation of the Hamiltonian is the average energy