Transcript Presentation Slides

ZCE 111 Computational Physics
Semester Project
Ng Wei Chun, 105525
Solving the Schrödinger equation for
Helium Atom with
Hatree-Fock Method
Helium Ground State with
Gaussian s-basis Function
Hamiltonian of a system of N
electrons, K nuclei, with Zn charges
• 𝐻=
1 1
4𝜋𝜀0 2
1 1
4𝜋𝜀0 2
𝑝𝑖2
𝑁
𝑖=1 2𝑚
+
𝑃𝑛2
𝐾
𝑛=1 2𝑀
𝑛
+
𝑒2
1
𝑁
𝑖,𝑗=1;𝑖≠𝑗 𝑟 −𝑟 − 4𝜋𝜀
0
𝑖
𝑗
2
𝑍
𝑍
𝑒
′
𝑛
𝑁
𝑛
𝑖,𝑗=1;𝑖≠𝑗 𝑅 −𝑅 ′
𝑛
𝐾
𝑛=1
𝑍𝑛 𝑒 2
𝑁
𝑖=1 𝑟 −𝑅
𝑛
𝑖
+
𝑛
(4.1)
• Index i for electrons, n for nuclei, m = electron
mass, Mn = nuclei masses.
Born-Oppenheimer approximation
• Nuclei are much heavier that they move much
slower than electron.
• 𝐻𝐵𝑂 =
−
1
4𝜋𝜀0
𝑝𝑖2
1 1
𝑁
𝑖=1 2𝑚 + 4𝜋𝜀 2
0
2
𝑍
𝑒
𝑛
𝑁
𝐾
𝑛=1 𝑖=1 𝑟 −𝑅
𝑛
𝑖
2
𝑒
𝑁
𝑖,𝑗=1;𝑖≠𝑗 𝑟 −𝑟
𝑖
𝑗
(4.2)
The independent-particle (IP)
Hamiltonian
• Further approximation reduce the 2nd term of
(4.2) to an uncoupled Hamiltonian.
• 𝐻𝐼𝑃 =
𝑝𝑖2
𝑁
𝑖=1 2𝑚
+ 𝑉 𝑟𝑖
(4.3)
Anti-symmetric trial wave function for
ground state
• Ψ 𝑥1 , 𝑥2 = Ψ 𝑟1 , 𝑠1 ; 𝑟2 , 𝑠2 =
1
𝜙(𝑟1 )𝜙(𝑟2 )
𝛼 𝑠1 𝛽 𝑠2 − 𝛼(𝑠2 )𝛽(𝑠1 )
2
(4.4)
• 𝛼 𝑠 = spin up, 𝛽 𝑠 = spin down, 𝜙 is an
orbital which share the same basis.
• The coordinates of the wave function are 𝑥1
and 𝑥2 , 𝑥𝑖 = 𝑟𝑖 , 𝑠𝑖 .
Born-Oppenheimer Hamiltonian for
the helium atom
• 𝐻𝐵𝑂 =
1
− 𝛻1 2
2
1
+ 𝛻2 2
2
+
1
𝑟1 −𝑟2
−
2
𝑟1
−
2
𝑟2
(4.5)
• Acting on (4.4),
•
1
2
− 𝛻1
2
+
1
2
𝛻2
2
−
2
𝑟1
−
2
𝑟2
+
• multiply both side by 𝜙 ∗ (𝑟1 ) and integrate over
𝑟2 :
•
1
2
− 𝛻1
2
−
2
𝑟1
3
+ 𝑑 𝑟2 𝜙 𝑟2
= 𝐸 ′ 𝜙(𝑟1 )
2
1
𝑟1 −𝑟2
𝜙 𝑟1
(4.7)
• All the integrals which yield constant multiple of
𝜙 𝑟1 are absorbed into 𝐸 ′ .
• Hereby, we have used 𝜙 ∗ 𝑟𝑖 𝜙 𝑟𝑗 = 𝛿𝑖𝑗 .
• Equation (4.7) is self-consistent.
Basis functions
• Beyond restricting the wave function to be
uncorrelated; by writing it as a linear
combination of 4 fixed, real basis functions.
• 𝜙 𝑟 =
4
𝑝=1 𝐶𝑝 𝜒𝑝
𝑟
(4.12)
• Leads directly from (4.7) that:
•
1
2
− 𝛻1 2 −
2
𝑟1
+
• Multiply 𝜒𝑝 𝑟1 from the left and integrate over 𝑟1 :
•
𝑝𝑞
′
𝐶
𝐶
𝑄
𝐶
=
𝐸
𝑟𝑠 𝑟 𝑠 𝑝𝑟𝑞𝑠 𝑞
ℎ𝑝𝑞 +
𝑝𝑞 𝑆𝑝𝑞 𝐶𝑞
(4.14)
• ℎ𝑝𝑞 =
1 2
𝜒𝑝 − 𝛻
2
• 𝑄𝑝𝑟𝑞𝑠 =
3
2
−
𝑟
𝜒𝑞
3
𝑑 𝑟1 𝑑 𝑟2 𝜒𝑝 𝑟1 𝜒𝑟 𝑟2
1
𝑟1 −𝑟2
𝜒𝑞 𝑟1 𝜒𝑠 𝑟2
• 𝑆𝑝𝑞 = 𝜒𝑝 𝜒𝑞 , overlap matrix
(4.15)
Gaussian s−basis function
• 𝜒𝑝 𝑟 = 𝑒
−𝛼𝑝 𝑟 2
(4.16)
• 𝑄𝑝𝑟𝑞𝑠 =
2𝜋5 2
𝛼𝑝 +𝛼𝑞 𝛼𝑟 +𝛼𝑠 𝛼𝑝 +𝛼𝑞 +𝛼𝑟 +𝛼𝑠
(4.17)
• 𝛼1 = 0.298073 ,𝛼2 = 1.242567 , 𝛼3
= 5.782948, and 𝛼4 = 38.474970
Program flow
• The 4 × 4 matrices ℎ𝑝𝑞 , 𝑆𝑝𝑞 and 4 × 4 × 4 × 4 array
𝑄𝑝𝑟𝑞𝑠 are calculated.
• Initial values for 𝐶𝑝 are chosen.
• Use C-values to construct matrix F:
• 𝐹𝑝𝑞 = ℎ𝑝𝑞 + 𝑟𝑠 𝑄𝑝𝑟𝑞𝑠 𝐶𝑟 𝐶𝑠
(4.18)
• Always check that 4𝑝,𝑞=1 𝐶𝑝 𝑆𝑝𝑞 𝐶𝑞 = 1, as the
normalisation.
(4.19)
• Solve for generalised eigenvalue problem
• 𝑭𝑪 = 𝐸 ′ 𝑺𝑪
•
•
•
•
(4.20)
For ground state, the vector C is the one corresponding to the
lowest eigenvalue.
Solution C from (4.20) is used to build matrix F again and so on.
Ground state energy can be found by evaluating the expectation
value of the Hamiltonian for the ground state:
𝐸𝐺 = 2 𝑝𝑞 𝐶𝑝 𝐶𝑞 ℎ𝑝𝑞 + 𝑝𝑞𝑟𝑠 𝑄𝑝𝑟𝑞𝑠 𝐶𝑝 𝐶𝑞 𝐶𝑟 𝐶𝑠
(4.21)
General Hartree-Fock Method
Particle-exchange operator, 𝑷𝑖𝑗
• 𝑷𝑖𝑗 Ψ 𝑥1 , 𝑥2 , … , 𝑥𝑖 , … 𝑥𝑗 , … , 𝑥𝑁
= Ψ 𝑥1 , 𝑥2 , … , 𝑥𝑗 , … 𝑥𝑖 , … , 𝑥𝑁
(4.22)
• For the case of an independent-particle
Hamiltonian, we can write the solution of the
Schrödinger equation as a product of oneelectron states:
• Ψ 𝑥1 , … , 𝑥𝑁 = 𝜓1 𝑥1 … 𝜓𝑁 𝑥𝑁
(4.24)
Spin-orbitals permuted
• The same state as (4.24) but with the spin-orbitals
permuted, is a solution too, as are linear combinations
of several such states,
• Ψ𝐴𝑆 𝑥1 , … , 𝑥𝑁
1
=
𝜖𝑃 𝑷𝜓1 𝑥1 … 𝜓𝑁 𝑥𝑁
𝑁! 𝑷
(4.26)
• 𝑷 is a permutation operator which permutes the
coordinates of the spin-orbitals only, and not their
labels, or acting on labels only (acting on one at a
time).
Slater determinant
• We can write (4.26) in the form of a Slater
determinant:
• Ψ𝐴𝑆 𝑥1 , … , 𝑥𝑁
𝜓1 𝑥1 𝜓2 𝑥1
⋯ 𝜓𝑁 𝑥1
1 𝜓1 𝑥2 𝜓2 𝑥2
⋯ 𝜓𝑁 𝑥𝑁
=
⋮
⋮
⋱
⋮
𝑁!
𝜓1 𝑥𝑁 𝜓2 𝑥𝑁 ⋯ 𝜓𝑁 𝑥𝑁
(4.27)
Hartree equation
• Fock extended the Hartree equation by taking anti-symmetry into
account.
• ℱ𝜓𝑘 = 𝜖𝑘 𝜓𝑘
(4.30)
• ℱ𝜓𝑘 = −
1
2
𝛻2
−
𝑁
𝑍𝑛
𝑛 𝑟−𝑅
𝑛
𝜓𝑘 𝑥 +
𝑑𝑥 ′ 𝜓𝑙 ∗
−
𝑙=1
𝑥′
𝑁
𝑙=1
𝑑𝑥 ′
𝜓𝑙
𝑥′
2
1
𝑟−𝑟 ′
𝜓𝑘 𝑥
1
′ 𝜓 𝑥
𝜓
𝑥
𝑘
𝑙
𝑟 − 𝑟′
(4.31)
• The operator ℱ is called the Fock operator.
Linear combinations of a finite
number of basis states 𝜒𝑝
• Expanding the spin-orbitals 𝜓𝑘 as linear
combinations of a finite number of basis states
𝜒𝑝 :
• 𝜓𝑘 𝑥 =
𝑀
𝑝=1 𝐶𝑝𝑘
𝜒𝑝 𝑥
(4.55)
• Then (4.30) assumes a matrix form
• 𝑭𝑪𝒌 = 𝝐𝒌 𝑺𝑪𝒌
(4.56)
• where 𝑺 is the overlap matrix for the basis used.
General form of the Fock operator
• ℱ =ℎ+𝐽−𝐾
(4.59)
• 𝐽 𝑥 𝜓 𝑥 =
• 𝐾 𝑥 𝜓 𝑥 =
∗
′
′
𝑑𝑥 𝜓𝑙 𝑥 𝜓𝑙 𝑥
𝑙
𝑙
′
∗
′
𝑑𝑥 𝜓𝑙 𝑥 𝜓 𝑥
′
1
𝜓
𝑟12
′
1
𝜓𝑙
𝑟12
𝑥
𝑥
(4.60)
• The sum over 𝑙 is over all occupied Fock levels.
• If written in terms of the orbital parts only, the Coulomb and
exchange operators read
• 𝐽 𝑟 𝜓 𝑟 =2
• 𝐾 𝑟 𝜓 𝑟 =
𝑙
𝑙
3 ′
∗
′
𝑑 𝑟 𝜙𝑙 𝑟 𝜙𝑙 𝑟
𝑑 3 𝑟 ′ 𝜙𝑙 ∗ 𝑟 ′ 𝜙 𝑟 ′
′
1
𝜙 𝑟
𝑟 ′ −𝑟
1
𝜙𝑙 𝑟
𝑟 ′ −𝑟
(4.61)
• In contrast with (4.60), the sums over 𝑙 run over half the
number of electrons because the spin degrees of freedom
have been summed over.
Fock operator
• The Fock operator now becomes
• ℱ 𝑟 = ℎ 𝑟 + 2𝐽 𝑟 − 𝐾 𝑟
(4.62)
• The corresponding expression for the energy is
found analogously and is given by
• 𝐸𝑔 = 2 𝑘 𝜙𝑘 ℎ 𝜙𝑘 + 𝑘(2 𝜙𝑘 𝐽 𝜙𝑘
Roothaan equation
• For a given basis 𝜒𝑝 𝑟 , we obtain the
following matrix equation, which is known as
the Roothaan equation:
• 𝑭𝑪𝑘 = 𝜖𝑺𝑪𝑘
(4.64)
• 𝑺 is the overlap matrix for the orbital basis
𝜒𝑝 𝑟
Matrix 𝑭
• 𝐹𝑝𝑞 = ℎ𝑝𝑞 +
𝑘
∗
𝐶
𝑟𝑠 𝑟𝑘 𝐶𝑠𝑘 2 𝑝𝑟 𝑔 𝑞𝑠 − 𝑝𝑟 𝑔 𝑠𝑞
(4.65)
• ℎ𝑝𝑞 = 𝑝 ℎ 𝑞 =
𝑑3 𝑟𝜒𝑝 ∗
𝑟
1 2
− 𝛻
2
−
𝑍𝑛
𝑛 𝑅 −𝑟
𝑛
𝜒𝑞 𝑟
(4.66)
•
𝑝𝑟 𝑔 𝑞𝑠 =
𝑑3 𝑟1 𝑑3 𝑟2 𝜒𝑝 ∗
∗
𝑟1 𝜒𝑟 𝑟2
1
𝑟1 −𝑟2
𝜒𝑞 𝑟1 𝜒𝑠 𝑟2
(4.67)
• 𝑘 labels the orbitals 𝜙𝑘 and p, q, and s label the basis
functions.
Density matrix
• The density matrix for restricted Hartree-Fock is
defined as
∗
• 𝑃𝑝𝑞 = 2 𝑘 𝐶𝑝𝑘 𝐶𝑞𝑘
(4.68)
• Using (4.68), the Fock matrix can be written as
• 𝐹𝑝𝑞 =
1
ℎ𝑝𝑞 +
2
𝑟𝑠 𝑃𝑠𝑟
2 𝑝𝑟 𝑔 𝑞𝑠 − 𝑝𝑟 𝑔 𝑠𝑞
(4.70)
Energy
• 𝐸=
1
+
2
𝑝𝑞 𝑃𝑝𝑞 ℎ𝑝𝑞
𝑃𝑝𝑞 𝑃𝑠𝑟
𝑝𝑞𝑟𝑠
1
𝑝𝑟 𝑔 𝑞𝑠 − 𝑝𝑟 𝑔 𝑠𝑞
2
(4.71)
Gaussian function
• 𝜒𝛼 𝑟 = 𝑃𝑀 𝑥, 𝑦, 𝑧 𝑒 −𝛼
𝑟−𝑅𝐴
2
(4.78)
• For 𝑙 = 0, these functions are spherically
symmetric, 1s-orbital is given as
• 𝜒𝛼 (𝑠) 𝑟 = 𝑒 −𝛼
𝑟−𝑅𝐴
2
(4.81)
The two-electron integral
•
1𝑠, 𝛼, 𝐴; 1𝑠, 𝛽, 𝐵 𝑔 1𝑠, 𝛾, 𝐶; 1𝑠, 𝛿, 𝐷
=
2
𝑑 3 𝑟1 𝑑 3 𝑟2 𝑒 −𝛼 𝑟1 −𝑅𝐴 𝑒 −𝛽 𝑟2−𝑅𝐵
2𝜋
=
2
1 −𝛾 𝑟 −𝑅 2 −𝛿 𝑟 −𝑅 2
1
𝐶 𝑒
2
𝐷
𝑒
𝑟12
5 2
𝛼+𝛾 𝛽+𝛿 𝛼+𝛽+𝛾+𝛿
𝛼𝛾
2
𝛽𝛿
× 𝑒𝑥𝑝 −
𝑅𝐴 − 𝑅𝐶 −
𝑅𝐵 − 𝑅𝐷
𝛼 + 𝛾𝛾
𝛽+𝛿
𝛼+𝛾 𝛽+𝛿
2
× 𝐹0
𝑅𝑃 − 𝑅𝑄
𝛼+𝛽+𝛾+𝛿
2
(4.123)
• 𝑅𝑃 =
𝑅𝐴 −𝑅𝐶
2
, 𝑅𝑄 =
𝑅𝐵 −𝑅𝐷
2
, 𝐹0 𝑥 = 𝑒𝑟𝑓 𝑥 =
′2
2 𝑥
′
−𝑥
𝑑𝑥 𝑒
𝜋 0