氫分子和水分子基態能量的Ab Initio計算

(1)

國

立

交

通

大

學

應用數學系

碩

士

論

文

氫分子和水分子基態能量的 Ab Initio 計算

Ab Initio computation for ground state energies of

.

H

2

and

H O

2

研究生：陳相宇

指導教授：葉立明教授

中

(2)

氫分子和水分子基態能量的 Ab Initio 計算

Ab Initio computation for ground state energies of

.

H

2

.

and

.

H O

2

研究生：陳相宇 Student：Shiang-Yen Chen

指導教授：葉立明 Advisor：Li-Ming Yeh

國立交通大學

應用數學系

碩士論文

A Thesis

Submitted to Department of Applied Mathematics

College of Science

National Chiao Tung University

in partial Fulfillment of the Requirements

for the Degree of

Master

in

Applied Mathematics

July 2009

Hsinchu, Taiwan, Republic of China

(3)

氫分子和水分子基態能量的 Ab Initio 計算

學生：陳相宇

指導教授

：

葉立明教授

國立交通大學應用數學學系﹙研究所﹚碩士班

摘

要

以 Ab initio 方法來計算氫分子和水分子的基態能量。我們計算能量的理論基

礎是密度泛函理論(DFT)和局部密度近似法(LDA)，我們首先以原子軌道為基底

函數，並對他們做正交歸一化，接著再求每一項能量密度泛函的值，最後運用

Lagrange multipliers 對總密度能量泛函做最小化，再以最速下降法求得基態總能。

(4)

Ab Initio computation for ground state energies

of

H

₂

and

H O

₂

student：

Shiang-Yen Chen

Advisors：Dr.

Li-Ming Yeh

Department of

Applied Mathematics

National Chiao Tung University

ABSTRACT

Computation of the ground-state energies of

H and

₂

H O is concerned. Based on the density

₂

functional theory (DFT) the ground state energy of a many-electron system can be expressed by

electronic density and is the minimum value of a density functional. So the computation can be done

by finding a minimum value of a total energy density functional. Firstly, we use linear combination

of atomic orbitals (LCAO) to construct the molecular orbital and calculate each term of the density

functional, for example, kinetic energy, electron-electron Coulomb energy, electron-core Coulomb

energy, and exchange-correlation energy. Next, we take an initial guess and minimize the total

density energy functional by using Largrange multipliers method. Finally by steepest descent method,

the ground state total energy is obtained. In some special care is needed for the computation of

Coulomb potential to avoid the singular points.

(5)

iii

誌

謝

誌

謝

誌

謝

誌

謝

首先我要感謝我的指導教授-葉立明老師。感謝老師每次都細心的為我解釋程

式上或是推導上不懂的地方，並且給我許多建議和解決問題的方法，甚至之後為

了讓我趕的上口試時間，還在回國的當天晚上特地回學校和我討論程式結果，並

且花了好多時間幫我修正我的論文，感謝老師煞費苦心的幫助我，讓我能順利完

成碩士學位。我還要感謝清大物理系-許貞雄老師，感謝許老師開了許多物理課

程，而且還在課後抽空替我解答問題，使我能對我要做的問題有更深的了解。

接著我要感謝我的爸爸媽媽，他們資助我讓我在生活上不用擔憂，感謝我的

哥哥在我感到無聊時陪我聊天，並且幫我修正我的論文中的英文語法。還要感謝

我的家人們，他們支持我，讓我在大學畢業後能繼續念研究所，而且沒有給我任

何壓力。

我還要感謝我的朋友們，感謝廖茜婷在我趕論文的時候陪我熬夜並且幫我修

正論文中的英文文法。感謝林浙于和楊長銘在這段時間裡，和我一起讀書、吃飯

和打電動、並且讓我接觸到基金和股票，使我對理財方面有更深的認識。感謝陳

廷彰在預官考試還有論文寫作上提供協助，沒有你們，我的研究所生活一定過的

相當乏味。

(6)

Abstract (in Chinese) i

Abstract (in English) ii

Acknowledgement iii

Contents iv

Figures vi

1 . P h y s i c s 1

1.1 Density functional theory (DFT)...1

1.2 Local density approximation (LDA)

………...3

1.3 Linear combination of atomic orbitals (LCAO)

...5

2 . P r o c e s s o f s t u d y

… … … …

. . . .

6

2.1 Orthonormalize basis functions………7

2.2 Kinetic energy functional

E

kin

[ ]

n

...8

2.3 .

_{Electron-electron Coulomb energy functional}

E

ee

[ ]

n

...9

2.4 Electron-nuclear Coulomb energy functional

E

ec

[ ]

n

...10

2.5 Exchange-correlation energy functional

E

xc

[ ]

n

...12

3. The iterative procedure

...1412

13

3.1 The gradient of the total energy functional

F

(C)

_...14

(7)

3.3 Numerical results of

H

2

………..18

3.4 Numerical results of

H Ｏ

2

_………26

3.5 Idea of pseudo potential………30

A Appendix

……….

A.1 The contents of the main program

……

….………31

A.2 The program for solving Poisson’s equation……..……86

Reference………166

(8)

Figures

Figure 1 The position of atoms in

H O

2

7 Figure 2 Extended interval N = 1 11

Figure 3 Extended interval N = 2 11

Figure 4

The procedure of computation 17

Figure 5

The graph of pseudo potential and pseudo wave function 30

(9)

1. _Physics

Nearly all physical properties are related to total energies or to the differences between

total energies. There is an advantage of Ab Initio total energy calculations. While just one

piece of given data is needed to calculate all the physical properties that are related to total

energy. The ground-state energy is meaningful in physics which

.

is

a stationary

.

state

.

with

the

lowest energy of a system and it is the most often state of a stable system.

.

The object of this

study is

to find

.

the

.

ground state total energies of

.

H

2

and

H O

2

by Ab Initio method.

.

...

1.1 Density-functional theory (DFT)

The first mentioned by Thomas and Fermi (1927), they submit a theory said we can use

electron density to express energy. But it is coarse of the era even there are some problems in

molecular system. The most famous people are Hohenberg and Kohn (1964), they submit the

Hohenberg - Kohn (H-K) theorem which contains two main statements of a many-electron

system. One is that the ground state properties of a many-electron system are uniquely

determined by an electronic density

n r

( )

. And the other is that a energy functional for the

system can be defined and it can be proved that the ground state electronic density minimizes

this energy functional. Therefore, we can put any guess density

n r

( )

into the total energy

functional

E[ ( )]

n r

until generate the minimum value which guarantees the ground-state

(10)

Unfortunately, Density-functional theory (DFT) only proves the existence of the

ground-state total energy functional

E [ ( )]

_G

n r

but it does not provide the exact form of the

total energy functional. In general, the ground-state total energy functional of a many-electron

system can be expressed by

E [ ( )] =

G

n r

E

ec

[ ( )]

n r

+

T n r

m

[ ( )]

+

E

H

[ ( )]

n r

where

E

ec

[ ( )]

n r

=

∫

v

ext

( ) ( )

r n r dr

and

v

ext

( )

r

is the external potential (the Coulomb

potential is due to the nuclear charges),

T [ ( )]

_m

n r

_{is the kinetic energy and}

E [ ( )]

_H

n r

is the

energy of electron-electron interaction. The last two terms of the density functional are

unknown.

One year later, Kohn-Sham (1965) provided a method to derive exact forms for

T [ ( )]

_m

n r

and

E [ ( )]

_H

n r

. They transform the intractable many-body problem of

interacting electrons in external potential to a non-interacting electrons in an effective

potential. Thus they can rewrite the total energy functional as the sum of the kinetic energy of

single-particle

T n r

₀

[ ( )]

, the Coulomb energy between electrons

E n r

_ee

[ ( )]

, and the rest

( called exchange-correlation energy

E

_xc

[ ( )]

n r

). Therefore, we have the form of the total

energy functional

………..………

G 0

(11)

1.2 Local density approximation (LDA)

The wave functions of electronic orbitals of atomic system are written as :

( )

( , )

i

r

R

pl

r Y

lm

φ

=

θ φ

represented by the index i={p,l,m}, the main quantum number p, the angular momentum

quantum number l, and

_{the magnetic quantum number m. The electron density is}

2

( )

_i _i

( )

i occupied

n r

f

φ

r

∈

=

∑

For example, O-atom has atomic number of Z =8 and occupies=1S,2S,2P, the

density of O-atom is

2 2 2 2 2 2

1S 1S 2S 2S 2P 2P 1S 2S 2P

( )

( ) +

( ) = 2

( ) +2

( ) +4

( )

n r

=

f

φ

r

f

φ

r

f

φ

r

φ

r

φ

r

φ

r

In this thesis, we use the atomic unit with

ℏ

=

m

= =

e

1 _.

Kohn and Sham (1965) provided a simplest method to treatment

E [ ( )]

_G

n r

. The total

energy functional is

G 0

E [ ( )] =

n r

E n r

_ee

[ ( )] [ ( )]

+

T n r

+

E n r

_ec

[ ( )]

+

E

_xc

[ ( )]

n r

1 ( ') ( )

=

'

2 '

n r n r

drdr

r

−

r

∫∫

+

T n r

₀

[ ( )]

+

E n r

ec

[ ( )]

+

E

xc

[ ( )]

n r

………

=

1 ₂

∫

v

ee

( ) ( )

r n r dr

+

T n r

0 [ ( )]

+

E n r

ec

[ ( )]

+

E

xc

[ ( )]

n r

(1.1)

Here, the electron density satisfies

N =

∫

n r dr

( ) ,

(1.2)

where N is the total electron number in this system.

E

_ee

[ ( )]

n r

is the electron-electron

Coulomb energy functional and

v

_ee

( )

r

is the electron-electron Coulomb potential

( )

ee

v

r

=

( ')

'.

'

n r

dr

r

−

r

∫

(12)

0 [ ( )]

T n r

is the kinetic energy functional expressed in terms of a system of

non-interacting electrons (under the assumption of independent particles)

………

2 * * 0

1 [ ( )]

( )

( ) =

( )

2

i i i i i i i occupied i occupied

T n r

f

φ

r

φ

r dr

f

φ

r

φ

r dr

∈ ∈



∇



=

_

−

_

∇

⋅∇





∑

∫

∑

∫

[ ( )]

ec

E

n r

is the electron-nuclear Coulomb energy functional

E

ec

[ ( )]

n r

=

∫

V

ext

( ) ( )

r n r dr

where

v

ext

( )

r

=

Z

r

−

_{is the external potential ( Z is an atomic number and the Coulomb}

potential is due to the nuclear charges) .

E

xc

([ ])

n

is the exchange-correlation energy

functional

E

xc

([ ])=

n

∫

ε

xc

( ( )) ( )

n r

n r dr

_,

where

ε

_xc

( ( ))

n r

is the exchange-correlation energy per unit density ( more detail of

( ( ))

xc

n r

ε

will be explained in section 2.5).

To compute molecular system, we should add one more term

E

cc

which is the core-

core Coulomb energy given by

,

1

2

I J cc I J I J

Z Z

E

τ

=

−

∑

.

(13)

1.3 Linear combination of atomic orbitals (LCAO)

Linear combination of atomic orbitals is a technique to calculate molecular orbitals.

A mathematical description is

1 1 2 2 n n

i i

,

i

C

ψ

=

φ

+

φ

+ …

φ

=

∑

φ

where

ψ

_{is a molecular orbital,}

φ

i

,

i

= …

1, , ,

n

are the atomic orbitals,

,

1, , ,

i

C

i

= …

n

are the coefficients corresponding to

φ

i

. Roughly speaking, we shall use

atom’s wave functions to construct the molecular wave function.

(14)

2. Process of study

In this study, we want to find the ground-state total energies of

H and

₂

H O . We take

₂

the atomic orbitals as the basis functions to construct the molecular orbital. Because these

basis functions are at different sites, we have to orthonormalize of them firstly. Next, we take

these orthonormal functions to calculate each term of the density functional, for example of

kinetic energy, electron-electron Coulomb energy, electron-core energy, and exchange-

correlation energy. Then, we write the total energy functional as

F

(C)

_where

1 2

C = ( , ,....,

C C

C

_n

)

are the coefficients which corresponded to the basis functions. Finally,

we take an initial guess and minimize the total energy functional

F

(C)

_{by using Lagrange}

(15)

2.1 Orthonormalize basis functions

We take the radial wave function

R

_pl

( )

r

of Hydrogen-like atom to be construct our

basis function. For example,

_{H has two H-atoms with distance 0.74Å. In Cartesian}

₂

coordinate,

.

we put the first H atom on the origin and the other H atom on point (0.74 , 0 , 0).

Each of them takes

.

five orbital

.

basis

.

functions which are written as

(1,0,0)

( )

r

φ

,

φ

(2,0,0)

( )

r

,

φ

(2,1,-1)

( )

r

,

φ

(2,1,0)

( )

r

,

φ

(2,1,1)

( )

r

, therefore, we use ten basis

functions to compute

H . In the same way,

2

H O has one O-atom and two H-atoms (the

2

position as illustrated in Figure 1.) O-atom takes six orbital basis functions and each of

H-atoms takes two orbital basis functions, therefore, we use ten basis functions to compute

2

H O

_.

Figure 1 : The position of atoms in

H O

₂

Because the atoms are at different

.

sites,

.

_{we use Gram-Schmidt method to these basis}

functions. By the Gram-Schmidt method,

………...………...

(16)

1 1

'( )

( )

i

( ) | ( ) ( )

i i j i j j

r

φ

− =

=

−

∑

〈

〉

.

'

''( )

'

i

r

φ

=

where

φ

i

'

=

1

2 ' | '

i

φ φ

〈

〉

_……...

For convenience, the orthonormalized basis functions are denoted by

1

( ) ( ), ,

r

2

r

10

( )

r

φ

,

φ

…

φ

.

Then, by linear combination of the above orthonormal basis functions, we can construct the

wave function by

( )

i i

( )

i

r

C

r

ψ

=

∑

φ

_(2.1)

The density of molecule is

2 2 _* ,

( )

( ) =

_{i i}

( ) =

_i _j

( ) ( )

_i _j i i j

n r

=

ψ

r

∑

C

φ

r

∑

C C

φ

r

φ

r

_(2.2)

2.2 Kinetic energy functional

E

kin

[ ( )]

n r

kin

E

₌

2

-( )

( )

2 r

r

ψ

∇

ψ

〈

|

〉

=

2 * ,

-( )

( )

2

i j i j i j

C C

〈

φ

r

|

∇

|

φ

r

〉

∑

=

* * ,

1 ( )

( )

2

i j i j i j

C C

∇

φ

r

⋅ ∇

φ

r dr

∑

_∫

.

………..………… .

Vector

r

=

( , , )

x y z

in rectangular coordinate

.

Then

* * * *

_{( , , ) (}

i

_,

i

_,

i

₎

i

x y z

x

y

z

φ

∂

∇

=

∂

(17)

( , , ) (

j

,

j

,

j

)

j

x y z

x

y

z

φ

∂

∇

=

∂

* * * *

_{( , , )}

i j i j i j i

x y z

j

x y z

x

y

z

φ

∂

∇

⋅ ∇

=

+

∂

……

We use the following space discretization :

* *

₍

_{, , )}

*

₍

_{, , )}

,

2

i i

x

h y z

i

x

h y z

x

h

φ

∂

+

−

=

∂

(

, , )

(

, , )

,

2

j j

x

h y z

j

x

h y z

x

h

φ

∂

+

−

=

∂

*

=

*

( ,

, )

*

( ,

, )

,

2

i i

x y

h z

i

x y

h z

y

h

φ

∂

+

−

∂

( ,

, )

( ,

, )

,

2

j j

x y

h z

j

x y

h z

y

h

φ

∂

+

−

=

∂

* *

_{( , ,}

₎

*

_{( , ,}

₎

,

2

i i

x y z

h

i

x y z

h

z

h

φ

∂

+

−

=

∂

( , ,

)

( , ,

)

.

2

j j

x y z

h

j

x y z

h

z

h

φ

∂

+

−

=

∂

The computation domain is

3

[a,b] ,

h

b

a

,

m

−

=

m

is the partition number on [a,b], (

m

must be even, since

.

we use

.

the Simpson’s rule ).

2.3 Electron-electron Coulomb energy functional

E

_ee

[ ]

n

1 ( ') ( )

[ ( )]

'

2 '

ee

n r n r

E n r

dr dr

r

=

−

∫∫

By (2.2),

* ,

( ') ( ')

( ')

_'

'

i j ee i j i j

r

n r

v

n r

dr

C C

dr

r

φ

( ( )) =

−

=

∑

−

∫

(18)

2 2

(

( ')) (

( ))

1 [ ( )]

'

2 '

i i k k i k ee

C

r

C

r

E

n r

dr dr

r

φ

=

−

∑

∫∫

* * , , ,

( ') ( ')

( ) ( )

'

i j k h i j k h i j k h

r

C C C C

dr dr

r

φ

=

−

∑

_∫∫

Now we need to calculate

*

_{( ') ( ') ( ) ( )}

*

1 '

( ') ( ')

( ) ( )

'

i j k h i j k h

r

dr dr

r

φ

=

−

∫∫

(2.3)

and that is a 6-dimension integral. By solving the Poisson’s equation :

( )

( ) ( )

( )

0

i j

U r

r

on

U r

on

φ

∆

= −

Ω





=

∂Ω



,

for

Ω

large enough, (2.3) is

∫

U r n r dr

( ) ( )

.

2.4 Electron-nuclear Coulomb energy functional

E n

ec

[ ]

E

ec

[ ( )]

n r

=

∫

v

ext

( ) ( )

r n r dr

In molecule,

( ) =

-I ext I I

Z

v

r

τ

−

∑

_{, where}

I

Z

is the I -th atom’s atomic number and

τ

I

is its atomic site.

2

[ ( )]

( )(

( ))

ec ext i i i

E

n r

=

∫

v

r

∑

C

φ

r

dr

* ,

_i _j _ext

( ) ( ) ( )

_i _j i j

C C

v

r

φ

r

φ

r dr

=

∑

∫

Thus

,

the

electron-nuclear Coulomb energy corresponding to each orbital can be written by

i

v

ext

r

j

φ

〈 |

( ) |

〉

_for

i j

_,

=

_1,2,

…,

n

(19)

some grid points may be close to

τ

I

and this results in computational inaccuracy.

We divide the integral range into two parts : A ; B. Part-A contains the singular point and

the rest is part-B. After partitioning the integral range

[a,b]

, finding the minimum interval

which contains the singular point named

[c ,d ].

0 0

_..

Next, we extend interval

[c ,d ] for N

0 0

intervals as [c ,d ]

_N _N

_{namely part-A. Then, taking a large partition number to part-A for}

getting a more precise result. Finally, we integrate part-A and part-B and sum them up.

For example, the extend interval N equals to 1 and 2 as illustrated in Figure 2 and Figure 3.

Figure 2 : Extended interval N = 1 Figure 3 : Extended interval N = 2

Then

〈 |

φ

i

v

ext

( ) | 〉

r

φ

j

can be expressed by

* * [ , ]\A-

A-( ) A-( ) A-( )

( ) ( ) ( )

ext i j ext i j a b part part

v

r

φ

r

φ

r dr

+

v

r

φ

r

φ

r dr

∫

(20)

2.5 Exchange-correlation energy functional

E n

xc

[ ]

E

_xc

([ ])

n

=

∫

ε

_xc

( ( )) ( )

n r

n r dr

* ,

( ( )) ( ) ( )

i j xc i j i j

C C

ε

n r

φ

r

φ

r dr

=

∑

∫

The exchange-correlation energy per unit density

ε

xc

( ( ))

n r

depends only on the

density

n r

( )

.

at the position

r

.

ε

xc

( )

n

can be decomposed into a sum of the exchange

energy

ε

x

( )

n

.

and the correlation energy

ε

c

( ).

n

( ) = ( )+ ( )

xc

n

x

n

c

n

ε

They are obtained from a system of uniform electron

.

gas with the corresponding

density

n

and they are frequently expressed in terms of the radius

r

s

of a unit charge

defined by

1 3

3 (

)

4

s

r

n

π

=

.

The exchange energy is

x

( )

0.458

s

n

r

ε

= −

and

the correlation energy is

( )

(1

1 2

) ;

1 ln

ln

;

1

s s c s s s s

r

for

n

A

r

B

Cr

r

Dr

for

γ

β

ε

= 

−



+



≧

＜

where

γ

= 0.1423 ;

β

₁

= 1.0529 ;

β

₂

= 0.3334

A = 0.311 ; B =

－

0.048 ; C =0.002 ; D =

－

0.0116 .

(21)

3. The iterative procedure

Because of the total energy density functional is

[ ( )]

tol

kin

ee

ec

xc

E

n r

=

E

n r

+

E

n r

+

E

n r

+

E

n r

and because the density

n r

( )

can be represented by

* ,

( )

_i

( ) ( )

_j _i _j i j

n r

=

∑

C C

φ

r

φ

r

_,

we can write the total energy functional as

F

(C)

_where

C ( , ...

=

C C

₁ ₂

C

₁₀

)

. Next, we

minimize the total energy functional

F

(C)

under the constraints :

* ,

N =

( )

_i _j _i _j i j

n r d r

=

∑

C C

φ φ

d r

∫

.

Because the basis functions

{ }

φ

_{i i}

_{= …}

_{1, ,10}

are orthonormal, the constraints can be

written as

C

12

+

C

22

+ … +

C

102

=

N. Then the problem becomes

1 2 10

min ( ,

F C C

…

C

)

under the constraints

C

₁2

+

C

₂2

+ … +

C

₁₀2

=

N

Define

G C C

( ,

₁ ₂

,…

,

C

₁₀

)

=

C

₁2

+

C

₂2

+ …+

C

₁₀2

N = 0

−

_{by Largrange multipliers,}

we need to solve the problem

2

1

2

10

1

2

10

2

1

2

10 N

or

(2 ,2

)

N

F

G

C

F

C

λ

∇ = ∇





+

+ … +

=



∇ =

…2





+

+ … +

=



Therefore, we have

(22)

1 1 2 2 10 10 2 2 2 1 2 10

2

2 N

F

C

F

C

F

C

λ

∂



=

∂



∂



=

∂





 ∂



=

∂



+

+ … +

=





⋮

(3.1)

.

3.1 The gradient of the total energy functional

F

(C)

The total energy functional

F C C

( ,

1 2

,…

,

C

10

)

can be rewrite by

(C)

_kin

(C)

_ec

(C)

_ee

(C)

_xc

(C)

F

=

E

+

E

+

E

+

E

_.

(3.2)

Thus

(

kin

(C)

ec

(C)

ee

(C)

xc

(C))

i i

F

E

C

∂

=

+

∂

.

From section 2.2~2.5, we can get

(

)

10 10 2 2 2 1 1

(C)

1 _{| -}

_|

_{| -}

_|

_{| -}

_|

2

kin j i j j i j i j j j i

E

C

=

φ

=

φ





∂

=

_

〈

∇

〉 + 〈

∇

〉

_

=

〈

∇

〉

∂

_

∑

_

∑

(3.3)

(

)

10 10 1 1

(C)

2

ec

j i ext j j ext i j i ext j

j j i

E

C

v

C

v

C

=

φ

=

φ





∂

=

_

〈 |

|

〉 + 〈 |

| 〉

_

=

〈 |

|

〉

∂

_

∑

_

∑

(3.4)

we note

1

1 '

'

i j k h k h i j

r

φ φ

= 〈

φ φ

〉

−

(23)

10 10 , , 1 , , 1 10 10 , , 1 , , 1

(C)

1

1 ₊

1

2 '

'

1

1 +

+

'

ee i j k i j k h i j h i j k h i j k i j h i i k h i j k h j k h i j k h i k h j k h

E

C C C

C

r

C C C

r

φ φ

= = = =



∂

=



_

∂

_

−





−

_

∑

10 10 , , 1 , , 1

1

1 '

'

i j k i j k h i j h i j k h i j k i j h

C C C

r

φφ

φ φ

φφ

φ φ

= =





=



+





₋





∑



10 , , 1

1

2 '

i j k i j k h i j k

C C C

r

φ φ

=

−

∑

(3.5)

If an initial guess

( ,

C C

₁ ₂

,…

,

C

₁₀

)

is given, then we can compute

〈 |

φ

i

v

ee

|

φ

j

〉

and (3.5)

becomes

10 1

2

_j _i _ee _j

.

j

C

φ

v

φ

=

〈 |

|

〉

∑

(3.6)

(

)

10 10 1 1

(C)

2

xc j i xc j j xc i j i xc j j j i

E

C

=

φ ε

φ

φ ε

φ

=

φ ε

φ





∂

=

_

〈 |

|

〉 + 〈 |

| 〉

_

=

〈 |

|

〉

∂

_

∑

_

∑

_(3.7)

To calculate (3.7), we need insert initial guess to decide

ε

xc

.

So we have

(

)

10 2 1

(C)

=

_j _i

-

_j

+2

_i _ext _j _i _ee _j _i _xc _j

.

j i

F

C

v

C

=

φ

φ ε

φ

∂

〈 | ∇ | 〉 〈 |

|

〉 + 2〈 |

|

〉 + 2〈 |

|

〉

∂

∑

Denote it by :

10 i j 1

(C)

=

j

A

j i

F

C

₌

∂

∑

(24)

10 1 j 1 1 10 2 j 2 1 10 10 j 10 1 2 2 2 1 2 10

A

2 A

2 N

j j j j j j

C

λ

= = =



=



=







₌



+

+ … +

=





∑

⋮

(3.8)

In matrix form (3.8) is

1 1 1 1 1 1 2 10 2 2 2 2 2 1 2 10 10 10 10 10 10 1 2 10 2 2 2 1 2 10

A

= 2

A

N

C

λ



 









 









 









 









_

 









_

_

_

_







₊

_{+ … +}

₌



⋯

⋮

⋱

⋮

⋯

(3.9)

3.2 The procedure of Ab Initio self-consistency

For a given

C

old

guess, we solve the eigenproblem (3.9), find the eigenvector

corresponding to the minimum eigenvalue, and let the eigenvector satisfy

2 2 2

1 2 10

N. C

+

C

+ … +

C

=

In this way, we obtain new coefficients

C

new

( , , ,

₁ ₂ ₁₀

)

C C

C

=

…

.

Compare new coefficients

C

new

with the previous one

C

old

C

< 1E-6 .

C

New Old New



₋















(25)

compute the total energy functional of coefficients

F

to get total energy (adding

E

cc

).

Otherwise we replace

C

old

by

C

new

as the initial guess

and repeat the above procedure , as

illustrated in Figure 4,

(26)

(27)

(28)

(29)

(

( experiment

experimental

experiment

al

al value of e

value of e

value of energ

nerg

nergyyyy for H

for H

2222

≒

≒ ----1.154948

≒

1.154948

1.154948 ～ ----1.175482 )

1.175482 )

The tables above do not include Exc. After observing the results given above and

comparing the total energies with each other, we can discover that the energy of Eec is the

main reason leading to an un-precise outcome. It shows that we could not get a more accurate

(30)

(31)

This time, we change the integral range of Ekin and Eee and also add Exc.

By observing the values of total energies of

.

result.1 and result.5, we can find that these total

energies do not change too much. By observing result.6, we can know that we get a even

more precise outcome, namely being more close to the experimental value, through doing

(32)

(33)

Final, we change the integral range of Eec and then observe the results. By comparing

with result.1, result.5 and result.7, we can discover that values of total energies have a big

change after changing the range of Eec. And also we can find that the total energies still do

(34)

(35)

(36)

(37)

((((experiment

experiment

experimental

al

al value of e

value of e

value of energ

nergyyyy for

nerg

for

H Ｏ≒

₂

≒

-75.58596)

Through

observing the results above, we can discover that they even do not converge.

From the results of

H and

₂

H Ｏ, we could know that it’s necessary to avoid the singular

2

points while calculating Eec, otherwise the outcome will not be precise even there is a

bigger partition number. It represents that we need to use other methods to calculate Eec. For

(38)

3.5 Idea of pseudo potential

The pseudo potential is used to represent the complicated effects of the movements of an

atom’s electrons and nucleus with other effective potentials. We should know Kohn-Sham

equation first. Kohn-Sham equation is obtained from the energy functional (1.1) by using the

variational principle with respect to

φ

*

( )

r

under the constraint

N =

∫

n r dr

( ) ,

2

1 ( )

( ( ))

( )

( ) ( )

2 v

ext

r

v

ee

n r

v

xc

n r

φ

i

r

ε

i

r

φ

i

r





− ∇ +

+

=









where

ε

i

is a Lagrange multiplier and eigenvalue corresponding to

φ

i

( )

r

.

Next, we need to construct a pseudo wave function which has the same graph as real

wave function while

r

>

r

c

₍

r

c

is a cutoff

.

radius) and is a smooth nodeless function to 0

as

r

< .

r

c

_{When the pseudo wave function have the same eigenvalue as real wave}

function, the potential is pseudo potential. Because the pseudo wave function is a smooth

nodeless function as

r

<

r

c

, there is no such kind of problems about singular point by using

pseudo potential.

(39)

A Appendix

A.1 The contents of the main program

!---!

! Ylm → 動能:complex 其他:real ! 積分 → 辛普森積分 !---宣告常數---!

module global implicit none

real(kind=8),parameter :: a0 = 0.529189379 ! Bohr radius 單位:Å real(kind=8),parameter :: pi = 3.141592653589 ! 圓周率π real(kind=8) :: x_max ! 積分範圍 x_min→x_max real(kind=8) :: y_max ! 積分範圍 y_min→y_max real(kind=8) :: z_max ! 積分範圍 z_min→z_max real(kind=8) :: x_min ! 積分範圍 x_min→x_max real(kind=8) :: y_min ! 積分範圍 y_min→y_max real(kind=8) :: z_min ! 積分範圍 z_min→z_max real(kind=8) :: EecX_min,EecY_min,EecZ_min real(kind=8) :: EecX_max,EecY_max,EecZ_max real(kind=8) :: Eec,Ecc,Eee,Eexc,E_kin,E_tol ! 各項能量 integer,parameter :: inn_m_parti=200 ! 內積積分-辛普森法切割數(必須為偶數) integer,parameter :: basics = 10 ! 基底函數數量 real(kind=8) :: C(basics) ! 基底函數係數 real(kind=8) :: C_new(basics),norm2,norm3 ! 基底函數係數 complex(kind=8) :: singleEkin(basics,basics) complex(kind=8) :: singleEec(basics,basics) complex(kind=8) :: singleEec2(6,basics,basics) complex(kind=8) :: ijkh_Eee(basics,basics,basics,basics) complex(kind=8) :: singleEee(basics,basics) complex(kind=8) :: singleExc(basics,basics) complex(kind=8) :: Exc_N real(kind=8) :: Exc_N2 real(kind=8)::xi(inn_m_parti+1),yi(inn_m_parti+1),zi(inn_m_parti+1) complex(kind=8),allocatable:: wave(:,:,:,:) ! Φi,i=1,10

(40)

real(kind=8),allocatable :: vext(:,:,:) real(kind=8),allocatable :: Exc(:,:,:)

integer,parameter :: CA_M =(1+BASICS)*BASICS/2 ! 計算上三角的總數 integer :: IC(CA_M),JC(CA_M) ! 計算上三角 end module !==============================主程式===============================! program main use IMSL use global implicit none !---宣告變數&常數---! real(kind=8) :: r,x1,y1,z1,x2,y2,z2 ! 座標變數 real(kind=8) :: theta,phi real(kind=8),external :: Rpl,vextF complex(kind=8),external :: Ylm complex(kind=8),external :: basis_func integer :: i,j,k,h,h2,counter,counter2 complex(kind=8) :: temA,temB,temC,temD

real(kind=8) :: F_diffC(basics,basics) ! F_diffC(basics)為 F 對 C 微分 real(kind=8) :: ANS

complex(kind=8) :: ANS2,ANS3,ANS4,ANS5 integer,parameter::nn=basics ! 計算 eigenvalue integer,parameter::nm=basics

integer is1,is2,ierr,matz

double precision a(nm,nn),wr(nn),wi(nn),z(nm,nn),fv1(nn),extremaN integer iv1(nn),select,extrema allocate(wave(basics,inn_m_parti+1,inn_m_parti+1,inn_m_parti+1)) allocate(grad_phi(inn_m_parti+1,inn_m_parti+1,inn_m_parti+1)) allocate(grad_phiA(inn_m_parti+3,inn_m_parti+3,inn_m_parti+3)) allocate(grad_phiB(inn_m_parti+3,inn_m_parti+3,inn_m_parti+3)) allocate(vext(inn_m_parti+1,inn_m_parti+1,inn_m_parti+1)) allocate(Exc(inn_m_parti+1,inn_m_parti+1,inn_m_parti+1)) !---把輸出結果印在 txt 檔上 open(unit=8,file='10wave_function.txt',status='old') ! 同檔案 10wave_function.txt open(unit=60,file='new_10single_Ekin.txt',status='old') ! 此檔案同 Ekin.txt

(41)

open(unit=61,file='new_10single_Eec.txt',status='old') ! 此檔案同 Eec.txt open(unit=62,file='new_10single_Eee.txt',status='old') ! 此檔案同 Eee.txt open(unit=90,file='10TOTAL_ENERGY.txt') ! 最後收斂後的 Total Energy !-前置作業 1---! !----FIRST START 要先算一次各基底函數對應的值 ---! !!! 給定積分範圍 Exc 跌代需要先給訂積分範圍 x_max = 10.0 ! 積分範圍 x_min→x_max y_max = 10.0 ! 積分範圍 y_min→y_max z_max = 10.0 ! 積分範圍 z_min→z_max x_min = -10.0 ! 積分範圍 x_min→x_max y_min = -10.0 ! 積分範圍 y_min→y_max z_min = -10.0 ! 積分範圍 z_min→z_max !----(0) 需用到正交歸一的波函數下列方法 2 選一 ---! !---(i) 對波函數正交歸一化 ! call othonormal() !---(ii)讀取先前正交歸一的波函數檔案 do h = 1,basics do i = 1,inn_m_parti+1 do j = 1,inn_m_parti+1 do k = 1,inn_m_parti+1 read(8,*) temD wave(h,i,j,k) = temD end do end do end do end do !----(1)動能部分 : 計算每個 <Φiα｜T｜Φjβ)> ---! ! call calEkin() ! !----(2)電-核庫倫能部分 : 計算每個 <Φiα｜Vext｜Φjβ)> ---! ! call calEec(2) !!給定奇異點區間擴展數 !----(3)電-電庫倫能部分 : 讀取 2 檔案存成<i,j| |k,h> ---! ! call calEee() ! GOTO 115 ! !---! !---- 讀檔案給定 singleEkin、singleEec 值 ---!

(42)

do i = 1,basics do j = 1,basics

read(60,'(TR8,(E30.20,E30.20))')temA ! TR8:從第 8 個位置開始讀 read(61,'(TR8,(E30.20,E30.20))')temB

singleEkin(i,j) = temA ! <Φj｜Ｔ｜Φi)> singleEec(i,j) = temB ! <Φj｜Vext｜Φi)> do k = 1,basics

do h = 1,basics

read(62,'(TR13,(E30.20,E30.20))')temC

ijkh_Eee(i,j,k,h) = temC ! <Φj｜Ｔ｜Φi)> end do end do end do end do ! 給定數值:(1)選取最小特徵值 (2) 選取最大特徵值 do select = 1,2 if(select==1)then write(90,*)"(1)選取最小特徵值:" write(90,*) elseif(select==2)then write(90,*) write(90,*)"(2)選取最大特徵值:" write(90,*) else write(*,*)"select 錯誤!" GOTO 115 endif !-- 給定初始係數 C(i)---! C(1) = 1.0 C(2) = 1.0 C(3) = 1.0 C(4) = 1.0 C(5) = 1.0 C(6) = 1.0

(43)

C(7) = 1.0 C(8) = 1.0 C(9) = 1.0 C(10)= 1.0 !---第一次迭代前的能量---! call calExc() ! 計算 Exc---

call calTol() ! 列出第一次的 total energy write(90,*) "第一次迭代前的能量!" write(90,*) "Ekin ", E_kin

write(90,*) "Eec " , Eec write(90,*) "Eee " , Eee write(90,*) "Exc " , Eexc write(90,*) "Total Energy ",E_tol write(90,*) counter = 1 ! 計算迭代係數 C 的次數 !--- 計算<Φi｜Vee｜Φj)> 存到 singleEee(:,:)，共 17*17 項 ---! 113 singleEee(:,:) = (0.0,0.0) ! 113 代表迭代係數從這裡開始 do i = 1,basics do j = 1,basics do k = 1,basics do h = 1,basics

! <Φi｜Vee｜Φj)> =ΣiΣkＣiＣj <Φk(r')Φh(r')｜Vee｜Φi(r)Φj(r))> singleEee(i,j) = singleEee(i,j) + C(k)*C(h)*ijkh_Eee(k,h,i,j) end do

end do end do end do

if(counter==1)GOTO 114

call calExc() ! 計算 <Φi｜εxc(n)｜Φj)>

!--- 給定 F_diffC(i,j) 實數 114 F_diffC(:,:) = 0.0

(44)

do j =1,basics

! F_diffC(i,j) = <Φi｜T + 2Vec + 2Vee + 2εxc｜Φj)>

F_diffC(i,j) = singleEkin(i,j)+2.0*singleEec(i,j) +2.0*singleEee(i,j)+2.0*singleExc(i,j) end do end do !----解 eigenvalue prob---! matz = 1 ! 給值 1 輸出特徵值&特徵向量 do i =1,basics do j =1,basics a(i,j) = F_diffC(i,j) end do end do call rg(nm,nn,a,wr,wi,matz,z,iv1,fv1,ierr) ! if (counter<=2)then ! 列出前 2 次跌代的特徵值與特徵向量 ! write(90,*) counter ! do j=1,nn ! write(90,*) ! write(90,*) "eigenvalne:",wr(j),wi(j) ! 列出特徵值實部虛部 ! write(90,*) "eigenvector:" ! do i=1,nn ! write(90,*) " ",z(i,j) ! 列出相對應的特徵向量 ! end do ! end do ! end if extremaN = wr(1) do j =1,nn if(select==1)then ! 取最小特徵值所對應的特徵向量 if(wr(j) <= extremaN )then

extremaN = wr(j) extrema = j end if

(45)

elseif(select==2)then ! 取最大特徵值所對應的特徵向量 if(wr(j) >= extremaN )then

extremaN = wr(j) extrema = j end if else write(*,*)"選取特徵值出錯!" end if end do ! !---取特徵值對應的特徵向量整理後令為新係數 C_new norm2 = 0.0 do i =1,basics

norm2 = norm2 + z(i,extrema)*z(i,extrema) end do do i =1,basics C_new(i) = (z(i,extrema)*sqrt(10.0))/sqrt(norm2) end do !---判停---! norm2 = 0.0 norm3 = 0.0 do i =1,basics

norm2 = norm2 + (C_new(i)-C(i))*(C_new(i)-C(i)) norm3 = norm3 + C_new(i)*C_new(i)

end do if(sqrt(norm2/norm3)<1E-6)then write(90,'(A2,I5,A8)') "第",counter,"次收斂!" C(:) = C_new(:) do i =1,basics write(90,*)" ", C(i) end do

call calTol() ! 計算 total energy else

(46)

write(90,*) counter,E_tol ! 列出每一次的 total energy do i =1,basics write(90,*)" ", C(i) end do counter = counter+1 C(:) = C_new(:) if(counter>100)then write(90,*) " 跌代次數超過 100 還不會收斂!" GOTO 115 else GOTO 113 end if end if !---- 計算個別的能量 ---! write(90,*)

write(90,*) "Ekin = " , E_kin write(90,*)

write(90,*) "Eec = " , Eec write(90,*)

write(90,*) "Eee = " , Eee write(90,*)

write(90,*) "Exc = " , Eexc write(90,*) !---- 計算 Ecc 項 ---! !---- 1 ＺI ＺJ ---! !---- Ecc ＝ —— ΣΣ ————— ---! !---- 2 I J ｜τI-τJ｜ ---! Ecc = 2.0*(8.0/0.9584) + 1.0/(0.9584*sqrt(2.0+2.0*cos(75.55*pi/180.0))) write(90,*) "Ecc = " , Ecc

write(90,*)

!--- 列出 Ekin + Eee + Eec + Ecc 總合---! write(90,*) "Total Energy =",E_kin + Eec + Eee + Eexc + Ecc

(47)

write(90,*) end do ! ← select 115 DEALLOCATE(wave) ! 釋放可變動陣列記憶體 DEALLOCATE(grad_phi) DEALLOCATE(grad_phiA) DEALLOCATE(grad_phiB) DEALLOCATE(vext,Exc) stop end !================= 主程式結尾 =============================! !--- radial function Rpl(r) ---! real(kind=8) function Rpl(atom,p,l,r)

use global implicit none integer :: p,l ! p 主量子數 ; l 角動量子數;atom 原子序 real(kind=8) :: r,atom if (p==1 .AND. l==0)then ! R10 Rpl = 2.0*((atom/a0)**(1.5))*exp(-(atom/a0)*r) else if (p==2 .AND. l==0)then ! R20

Rpl = ((atom/(2.0*a0))**(1.5))*(2.0-(atom/a0)*r)*exp(-(atom/(2.0*a0))*r) else if (p==2 .AND.l==1)then ! R21

Rpl = (1.0/sqrt(3.0))*((atom/(2.0*a0))**(1.5))*((atom/a0)*r)*exp(-(atom/(2.0*a0))*r) else if (p==3 .AND. l==0)then ! R30

Rpl = (2.0/27.0)*((atom/(3.0*a0))**(1.5))*(27.0-18.0*((atom/a0)*r)+ & & (2.0*atom*atom/(a0*a0))*r*r)*exp(-(atom/(3.0*a0))*r) else write(*,*)" Rpl 超過 R30 !" end if return end

(48)

!--- 為計算上方便省略 phi 部分 → 實函數 ---! complex(kind=8) function Ylm(l,m,theta,phi)

use global implicit none

real(kind=8) :: theta,phi ! θ,φ

integer :: l,m ! l 角動量子數 ; m 磁量子數 if (l==0 .AND. m==0)then ! Y00

Ylm = 1.0/(sqrt(4.0*pi)) else if(l==1 .AND. m==-1)then ! Y1-1

Ylm = (sqrt(3.0/(8.0*pi)))*sin(theta)*(cos(phi)-csqrt((-1.0,0.0))*sin(phi)) else if(l==1 .AND. m==0)then ! Y10

Ylm = (sqrt(3.0/(4.0*pi)))*cos(theta) else if(l==1 .AND. m==1)then ! Y11

Ylm = -(sqrt(3.0/(8.0*pi)))*sin(theta)*(cos(phi)+csqrt((-1.0,0.0))*sin(phi)) else if(l==2 .AND. m==-2)then ! Y2-2

Ylm = (sqrt(15.0/(32.0*pi)))*((sin(theta))**(2.0))*exp(-2.0*csqrt((-1.0,0.0))*phi) else if(l==2 .AND. m==-1)then ! Y2-1

Ylm = (sqrt(15.0/(8.0*pi)))*(sin(theta))*(cos(theta))*exp(-csqrt((-1.0,0.0))*phi) else if(l==2 .AND. m==0)then ! Y20

Ylm = (sqrt(5.0/(16.0*pi)))*(3.0*(cos(theta))**(2.0)-1.0) else if(l==2 .AND. m==1)then ! Y21

Ylm = -(sqrt(15.0/(8.0*pi)))*(sin(theta))*(cos(theta))*exp(csqrt((-1.0,0.0))*phi) else if(l==2 .AND. m==2)then ! Y22

Ylm = (sqrt(15.0/(32.0*pi)))*((sin(theta))**(2.0))*exp(2.0*csqrt((-1.0,0.0))*phi) end if

return end

!--- 原始基底函數 Φ1～Φ10 ---( Φ = Rpl*Ylm )---! complex(kind=8) function basis_func(i,x1,y1,z1)

use global implicit none

integer :: i

real(kind=8) :: x1,y1,z1,atom ! x1,y1,z1 為直角坐標上變數

real(kind=8) :: r,r1,r2,r3 ! r1 中心為 H1 r2 中心為第一個 H2 real(kind=8) :: theta,theta1,theta2,theta3 ! 角度 θ:theta

(49)

real(kind=8) :: phi,phi1,phi2,phi3 ! 角度 φ:phi real(kind=8),external :: Rpl ! 函數 Rpl

complex(kind=8),external ::Ylm ! 函數 Ylm

r1 = sqrt(x1**2.0+y1**2.0+z1**2.0) ! r 換成直角坐標 x,y,z 中心為 O r2 = sqrt((x1-0.9584)**2.0+y1**2.0+z1**2.0) ! 中心為 H1

r3 = sqrt((x1+0.9584*cos(75.55*pi/180.0))**2.0+(y1-0.9584*sin(75.55*pi/180.0))**2.0+z1**2.0) ! 中心為H2

theta1 = atan2(sqrt(x1**2.0+y1**2.0),z1) ! ← arctangent(a/b) theta2 = atan2(sqrt((x1-0.9584)**2.0+y1**2.0),z1) theta3 = atan2(sqrt((x1+0.9584*cos(75.55*pi/180.0))**2.0+(y1-0.9584*sin(75.55*pi/180.0))**2.0),z1) phi1 = atan2(y1,x1) phi2 = atan2(y1,x1-0.9584) phi3 = atan2(y1-0.9584*sin(75.55*pi/180.0),x1+0.9584*cos(75.55*pi/180.0)) if(i<=6)then ! 波函數共 10 項，前 6 項中心 O atom = 8d0 r = r1 theta = theta1 phi = phi1 elseif(i==7.OR.i==8)then ! 第 7,8 項中心 H-1 atom = 1d0 r = r2 theta = theta2 phi = phi2 elseif(i==9.OR.i==10)then ! 第 9,10 項中心 H-2 atom = 1d0 r = r3 theta = theta3 phi = phi3 else write(*,*)"wave function 出問題" end if if(i==1)then

氫分子和水分子基態能量的Ab Initio計算

國

立

交

通

大

學

應用數學系

碩

碩

碩

碩

士

士

士

士

論

論

論

論

文

文

文

文