原子基態能量的計算

國

立

交

通

大

學

應用數學系

碩

士

論

文

原子基態能量的計算

Computation of the ground state energy of atoms

研 究 生：陳廷彰

指導教授：葉立明 教授

原子基態能量的計算

Computation of the ground state

energy of atoms

研 究 生：陳廷彰 Student：Ting-Jang Chen

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

國 立 交 通 大 學

應 用 數 學 系

碩 士 論 文

A Thesis

Submitted to Department of Computer and 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

中華民國九十八年七月

原子基態能的量計算

學生：陳廷彰

指導教授

：

葉立明 教授

國立交通大學

應用數學學系﹙研究所﹚碩士班

摘

要

基態能量在量子力學裡是代表最低能量。我們利用從 Kohn-Sham 方程式以及

local density approximation 中化簡出來之電子密度

ρ

的基態能量泛函

E_G[ ]ρ

來計算

原子的基態能量。如果要計算

E_G[ ]ρ

，我們必須先算出原子基態的波函數以及電

子密度，而電子密度是跟波函數有關的函數。所以第一步，要利用解 Kohn-Sham

方程式

HΨ = Ψε

找出波函數。Kohn-Sham 方程式是一個二階偏微分方程，在 Koh-

n-Sham 方程式中，

ε

的最小值所對應的函數

Ψ

，就是原子基態能量的波函數。為

了在計算上的方便，我們將

H

離散化，讓問題變成解特徵值的計算之後，利用自

洽來解出波函數。此篇計算所得到的數據和實際值的數據誤差不超過 15%。

中 華 民 國 九 十 八 年 七 月

i

Computation of the ground state

energy of atoms

Student: Ting-Jang Chen

Advisor: Dr. Li-Ming Yeh

Department of Applied Mathematics

National Chiao Tung University

Hsinchu, Taiwan, R.O.C.

Abstract

The ground state is the lowest-energy state of the quantum me-chanical system. We compute the ground state energy by using the ground state energy functional EG[ρ] of electronic density ρ, which is

deduced from the Kohn-Sham total-energy functional and local den-sity approximation. To compute the ground state energy of atoms, we have to compute each electronic density of the ground state of atoms. The electronic density of atoms is a function about the wave function of atoms, so the first step is to determine the wave func-tions of ground state of atoms by solving the Kohn-Sham equation

HΨ = εΨ. The Kohn-Sham equation is a problem of the second

or-der partial differential equation. The wave function of the ground state of atoms is the function Ψ corresponding to the minimal ε of the Kohn-Sham equation. For the convenience of solving the Kohn-Sham equation, we discretize the Hamiltonian H of the Kohn-Sham equation and the problem become an eigenvalue problem. In this computation, we determine the wave function of the ground state of atoms by self-consistency. The errors between the computation and the realistic values are less than 15%.

iii

誌

謝

能完成這篇論文，首先要感謝我的指導教授葉立明教授。

有耐心的引導我學習，糾正我錯誤的學習態度以及習慣，並

常常為了我的論文和我討論到深夜，真的是由衷的感謝您。

再來要感謝和我一起畢業的同學以及學弟妹，無條件的幫

我檢查論文，討論論文的內容，讓我的論文更加完整。

還要感謝我的女友劉翊婷，在我最困難最煩惱的時候總是

陪著我幫我解決問題，半夜陪著我準備口試，不停的幫我加

油打氣。

也感謝我的口試委員葉立明教授、王夏聲教授、黃聰明教

授、李華輪教授。

最後要感謝我的父母，時時掛念我的健康，感謝以上各位

的幫助以及支持，使我度過充實的研究生活。

Contents

Abstract (in Chinese) i

Abstract (in English) ii

Acknowledgement iii

Contents iv

1 Introduction 1

1.1 Density functional theory . . . 1

1.1.1 Electronic density . . . 1

1.1.2 Density functional theory and the ground state energy functional . . . 1

1.1.3 Local density approximation . . . 2

1.2 The ground state energy functional . . . 3

1.3 Atomic System . . . 6

2 Calculate the ground state energy of atoms 8 2.1 The process of the computation . . . 8

2.2 The integration and the differentiation . . . 8

2.3 Discretization of −1 2 d2_P nl(r) dr2 . . . 9 2.4 Discretization of l(l + 1) 2r2 i . . . 10

2.5 Discretization of effective potential . . . 11

2.5.1 Discretization of Coulomb potential . . . 11

2.5.2 Discretization of exchange-correlation potential . . . . 11

2.5.3 Discretization of external potential . . . 12

2.6 Normalization of wave function . . . 12

2.7 Self-consistency . . . 13

3 Results 16

A Reference 20

B Appendix 21

1

Introduction

1.1

Density functional theory

1.1.1 Electronic density

The starting point of this computation is the observation of Hohenbergy and Kohn (1964) [1] that electronic density ρ contains in principle all the infor-mation contained in a many-electron wave function. The electronic density

ρ of a many-electron system at point ~r is defined to be (In this paper, we use

the atomic unit with ¯h = me = e = 1)

ρ(r) = ∑

i∈occupied

fi|Ψi|2,

where fi are the occupation numbers of the orbitals denoted by i, and Ψi

is the wave function. It is convenient to employ the spherical coordinates

r, θ, φ, the wave function can be written as :

Ψi = Rnl(i)(r)Ylm(i)(θ, φ). (1.1)

The electronic orbitals are represented by the index i ={n, l, m}. n, l, m are the main quantum number, the angular momentum quantum number and the magnetic quantum number respectively.

1.1.2 Density functional theory and the ground state energy

func-tional

Density functional theory is a theory of quantum mechanics. With this the-ory, the properties of a many-electron system can be determined by using the functionals which is dependent on the electronic density. Hence the name density functional theory comes from the use of functionals of the electronic density. Hohenberg and Kohn observed [1] that the ground-state energy is a

functional of the electronic density ρ and can written as

EG= EG[ρ] = T0[ρ] + Eee[ρ] + Eex[ρ] + Exc[ρ],

where T0[ρ] is the kinetic energy functional, Eee[ρ] is the functional the

clas-sical Coulomb repulsion energy, Eex[ρ] is the energy of the electrons in the

external field of the nuclei, and Exc[ρ] is the exchange-correlation energy.

1.1.3 Local density approximation

The principle of local density approximation is to compute the exchange-correlation energies per particle. Hohenberg and Kohn provided some moti-vation for using approximate methods to describe the exchange-correlation energies as a function of a electronic density. The simplest method of de-scribing exchange-correlation energy is local density approximation. In local density approximation, the exchange-correlation energy of an electronic sys-tem is written as [7, 10] Exc = ∫ R3 ∫ εxc· ρ(r)drdω, withεxc= εx+ εc. Notation:

Here we define the integral

∫ R3 ∫ g(r)drdω for a function g(r) to be ∫ R3 ∫ g(r)drdω = ∫ _π 0 ∫ _2π 0 ∫ _∞ 0 g(r)r2sinφdrdθdφ = 4π ∫ _∞ 0 g(r)r2dr. (1.2) 2

1.2

The ground state energy functional

We use the ground state energy functional to compute the ground state energy of atoms. In local density approximation [10] the ground state energy of a many-electron system are expressed as a functional EG[ρ] of the electronic

density ρ(r). The ground state energy functional can be written as

EG[ρ] = F [ρ(r)] + Eex[ρ]. And Eex[ρ] = ∫ R3 ∫ vex· ρ(r)drdω = ∫ R3 ∫ −z rρ(r)drdω =−4π ∫ _∞ 0 z· r · ρ(r)dr, (1.3) where z is the atomic number and vex(r) is the external potential

vex(r) =−

z r.

The functional F [ρ(r)] is given by [10]

F [ρ(r)] = Eee[ρ] + T0[ρ] + Exc[ρ] = 1 2 ∫ R3 ∫ Vee(r)ρ(r)drdω + T0[ρ] + Exc[ρ] = 1 2 ∫ R3 ∫ ∫ R3 ∫ _ρ(r)ρ(r0) |r − r0_| dr0dωdrdω + T0[ρ] + Exc[ρ].

And we obtain that

EG[ρ] = ∫ R3 ∫ Vee(r)ρ(r)drdω + T0[ρ] + Exc[ρ] + ∫ R3 ∫ vex(r)ρ(r)drdω. (1.4)

The first term of the right hand side of (1.4) is the Coulomb energy of the electron charge densities, with Vee(r) being the electron-electron Coulomb

potential [1]: Vee(r) = ∫ R3 ∫ ρ(r0) 1 r >dr 0_{dω = 4π}∫ ∞ 0 ρ(r0) r 02 r >dr 0

(here we define the variable ”r >” be max(r, r0)). So the first term of the right hand side of (1.4) is

1 2 ∫ R3 ∫ (∫ R3 ∫ _ρ(r)ρ(r0₎ |r − r0_| dr0dω ) drdω = 1 2·(4π) 2∫ ∞ 0 ∫ _∞ 0 r2_r02 r > ρ(r)ρ(r 0_)dr0_dr. (1.5) 3

Exc[ρ] is the exchange-correlation energy. In the local density approximation, Exc[ρ] is written as Exc[ρ] = ∫ Ω ∫ εxc[ρ(r)]ρ(r)drdω = 4π ∫ _∞ 0 εxc[ρ(r)]· ρ(r) · r2dr. (1.6)

Now, we discuss the exchange-correlation energy εxcin (1.6). The

exchange-correlation energy may be decomposed into a sum of the exchange energy and the correlation energy [10].

εxc= εx+ εc.

εxc is expressed in terms of the radius rs of a unit charge defined by

rs = (

3 4πρ)

1 3.

The exchange term can be calculated exactly for a uniform electron gas, and is given by εx[ρ] =− 3 4πrs (9 4) 1 3 =−0.458 rs .

The correlation energy is given by

εc =      − γ 1 + β1√rs+ β2rs ; f or rs ≥ 1, low density A ln rs+ B + Crsln rs+ Drs ; f or rs < 1, high density ,

where the values of the constants are given by

γ = 0.1423; β1 = 1.0529; β2 = 0.3334,

and

A = 0.0311; B =−0.0480; C = 0.0020; D = −0.0116.

The term T0[ρ] is the kinetic energy functional. The kinetic energy

func-tional is expressed in terms of a system of noninteracting electrons. (Here the wave function Ψ(r, θ, φ) = Rnl(r)· Ylm(θ, φ). )

T0[ρ] = ∑ i∈occupied fi ∫ R3 ∫ Ψ∗_i(−∇ 2 2 )Ψidrdω. (1.7) 4

The natation that∇2 in spherical coordinate can be written as [4] ∇2 ₌ 1 r2 ∂ ∂r(r 2 ∂ ∂r) + 1 r2_sin2_θ ∂2 ∂φ2 + 1 r2_sinθ ∂ ∂θ(sinθ ∂ ∂θ).

And we obtain that

−1 2∇ 2_Ψ i = − 1 2 { 1 r2 ∂ ∂r(r 2 ∂ ∂rΨi) + 1 r2 [ 1 sinθ ∂ ∂θ(sinθ ∂Ψi ∂θ ) + 1 sin2_θ ∂2_Ψ i ∂φ2 ]} = −1 2 { 1 r2 ∂ ∂r(r 2 ∂ ∂rRnl(i))Ylm(i)+ 1 r2Rnl(i) [ 1 sinθ ∂ ∂θ(sinθ ∂Ylm(i) ∂θ ) + 1 sin2_θ ∂2_Y lm(i) ∂φ2 ]} . (1.8) Because of the identity [4]

1 sinθ ∂ ∂θsinθ ∂Ylm(i) ∂θ + 1 sin2_θ ∂2Ylm(i) ∂φ2 =−l(l + 1)Ylm(i),

we obtain the equation

−1 2∇ 2_Ψ i = − 1 2 { 1 r2 ∂ ∂r(r 2 ∂ ∂rRnl(i))Ylm(i)− 1 r2Rnl(i) [₁ r2 · l(l + 1) · Ylm(i) ]} = { − 1 2r2 d dr(r 2dRnl(i) dr )+ l(l + 1) 2r2 Rnl(i) } Ylm(i). (1.9) By (1.9) we have ∫ R3 ∫ Ψ∗_i(−1 2∇ 2 )Ψidrdω = ∫ R3 ∫ Ylm(i)Rnl(i) [ − 1 2r2 d dr(r 2dRnl(i) dr )+ l(l + 1) 2r2 Rnl(i) ] Ylm(i)drdω = ∫ R3 ∫ |Ylm(i)|2Rnl(i) [ − 1 2r2 d dr(r 2dRnl(i) dr )+ l(l + 1) 2r2 Rnl(i) ] drdω. (1.10) By (1.7), (1.10), and the normalization of Ylm (in section 2.6 we will

discuss the normalization) [7]

∫ π 0 ∫ 2π 0 |Ylm(i)| 2 sinφdθdφ = 1, 5

we get the equation T0 = ∑ i∈occupied fi {∫ _∞ 0 r2Rnl(i) [ − 1 2r2 d dr(r 2dRnl(i) dr ) + l(l + 1) 2r2 Rnl(i) ] dr } = ∑ i∈occupied fi {∫ _∞ 0 −1 2 Rnl(i)· r · (2 · dRnl(i) dr + d2_R nl(i) r2 )dr + ∫ _∞ 0 Rnl(i) l(l + 1) 2r2 Rnl(i)r 2 dr } . (1.11) By (1.2), (1.4), (1.5), (1.6), (1.11), we get the formula of ground state energy of atoms. EG[ρ] = 1 2 ∫ _∞ 0 ∫ _∞ 0 r2r02 r >ρ(r)ρ(r 0_{)· (4π)}2 drdr0 + ∑ i∈occupied fi {∫ _∞ 0 −1 2 Rnl(i)· r · (2 · dRnl(i) dr + d2_R nl(i) r2 )dr + ∫ _∞ 0 Rnl(i) l(l + 1) 2r2 Rnl(i)r 2_dr } + 4π ∫ _∞ 0 εxc[ρ(r)]· ρ(r) · r2dr −4π∫ ∞ 0 z· r · ρ(r)dr. (1.12)

1.3

Atomic System

In order to calculate the ground state energy functional. we must find each wave function of the atoms by solving the Kohn-Sham equation. The Kohn-Sham equation can be written as [1]:

( −1 2∇ 2 + vef f[ρ] ) Ψi(r) = εiΨi(r), (1.13)

where εi is a Lagrange multiplier and can be interpreted as the orbital energy

of the state represented by Ψi.

By (1.9) and (1.13), we have the equation

( −1 2 1 r2 d dr(r 2dRnl dr ) + l(l + 1) 2r2 Rnl ) Ylm+ vef fRnlYlm = enlRnlYlm. (1.14)

Multiplying both side of (1.14) by _Y1

lm, we obtain the equation

−1 2 1 r2 d dr(r 2dRnl dr ) + l(l + 1) 2r2 Rnl+ vef fRnl= enlRnl. (1.15) 6

The effective potential can be written as

vef f = vex[ρ(r)] + Vee[ρ(r)] + vxc[ρ(r)]. (1.16)

In section 1.2 we have known that

vex =− z r, Vee= ∫ R3 ∫ ρ(r0) 1 r >dr 0_.

The exchange-correlation potential vxc is given by [1]

vxc[ρ(r)] =

δ (εxc[ρ(r)]· ρ(r))

δρ(r) = εxc[ρ(r)] + ρ(r) dεxc

dρ ,

where εxc= εx+εxis written in section 1.2. For the convenience of computing

(1.15), we presume a function Pnl such that

Pnl(r) = rRnl(r).

Multiplying both sides of (1.15) by r, this equation may also be written as −1 2 d2Pnl(r) dr2 + l(l + 1) 2r2 Pnl(r) + vef f[ρ]Pnl(r) = enlPnl(r). (1.17) So ₍ −1 2 d2 dr2 + l(l + 1) 2r2 + vef f[ρ] ) Pnl(r) = enlPnl(r). (1.18)

The electronic density ρ(r) is written as

ρ(r) = ∑ i∈occupied fi|Ψi(r, θ, φ)|2 = ∑ i∈occupied fi   ∑l m=−l |Rnl(i)(r)· Ylm(i)(θ, φ)|2   = ∑ i∈occupied fi   ∑l m=−l |Rnl(i)(r)|2· |Ylm(i)(θ, φ)|2   = ∑ i∈occupied fi  _|Rnl(i)(r)|2_· l ∑ m=−l |Ylm(i)(θ, φ)|2  _.

For a closed shell system we have the identity that [2]

l

∑

m=−l

|Ylm(θ, φ)|2 = 1,

we obtain the equation

ρ(r) = ∑

i∈occupied

fi|Rnl(i)(r)|2.

In section 2, we will describe the discretization of the equation (1.18).

2

Calculate the ground state energy of atoms

2.1

The process of the computation

In the beginning of this computation, we give an initial guess of the elec-tronic density to determine the correct wave function and elecelec-tronic density by self-consistency (In section 2.7 we will introduce self-consistency). After finding the correct wave function and electronic density, we can calculate the ground state energy by using the ground state energy functional EG. We

choose the electronic density builded by the wave function of hydrogen [2] to be the initial guess of the computation of helium. After the computation of helium being finished, we can use the electronic density of helium to be the initial guess of the computation of lithium. Such like the process of the computation of helium, we can compute the ground state energy of the atoms one by one by using the preceding atom to be the initial guess.

2.2

The integration and the differentiation

By the functional (1.12), we know that if we want to the compute ground state energy, we must find the wave functions and electronic density. To find the wave functions, we will solve the equation (1.18). Now, we discuss the computational method and the discretization of this computation. In this computation, we try to compute in the numerical method. We take the

computational domain r = [0˚A, 10˚A] and divide it into 400 equal parts.

The domain is represented as

~ r =          r1 r2 .. . r399          , where ri = h· i, and h is 10˚A over 400, h = 10˚A 400.

Now, we define the integration and the differentiation in this computa-tion. First, we use the ”trapezoidal rule” to compute the integration in this computation. The domain has been discretized as the vector ~r.

So we defined the integration as

∫ _rk ri f (r)dr = f (ri) + f (ri+1) 2 · h + f (ri+1) + f (ri+2) 2 · h +... +f (rk−1) + f (rk) 2 · h = h· (f (ri) 2 + k_∑−1 j=i+1 f (rj) + f (rk) 2 ) , f or 0≤ i < k ≤ 399 (2.19) ∫ _∞ 0 f (r)dr≈ ∫ 10˚A 0˚A f (r)dr = h· (f (r1) 2 + 398 ∑ j=2 f (rj) + f (r399) 2 ). (2.20) And the differentiation in this computation is defined as

f0(ri) = df (x) dx ¯¯ ¯¯ ¯ x=ri = f (ri+1)− f(ri−1) ri+1− ri−1 , f or 0≤ i < k ≤ 399. (2.21) The differentiation of f (r1),f (r399) is defined as

f0(r1) = f (r2)− f(r1) r2− r1 , f0(r399) = f (r399)− f(r398) r399− r398 . (2.22)

2.3

Discretization of

−

1

2

d

2

P

nl

(r)

dr

2 The discretization of −1₂d2Pnl(r) dr2 is an approximation. 9

By the Taylor series P (x + h) = P (x) + P0(x)h + P 00 2!h 2 ₊P000 3! h 3_{+· · · ,} and P (x− h) = P (x) + P (x)(−h) + P 00 2! (−h) 2₊ P000 3! (−h) 3_{+· · · .} So that P (x + h) + P (x− h) = 2P (x) + P00h2+ O(h4). We have −P (x + h) + 2P (x) − P (x − h) h2 ≈ −P 00_, and we obtain −1 2P 00 _{≈ M} p =          1 h2 _2h−12 · · · 0 −1 2h2 1 h2 . .. ... .. . . .. ... _2h−12 0 · · · _2h−12 1 h2          399×399 . (2.23)

2.4

Discretization of

l(l + 1)

2r

2_i l(l+1) 2r2 i is discretized as Ml =           l(l+1) 2r2 1 0 l(l+1) 2r2 2 . .. 0 l(l+1)_2r2 399           399×399 . (2.24) 10

2.5

Discretization of effective potential

By (1.16) we have known that effective potential contains Coulomb potential, exchange-correlation potential and external potential. Now we discuss their discretization.

2.5.1 Discretization of Coulomb potential

The Coulomb potential

vee = ∫ R3 ∫ _ρ(r0) |r − r0_|dr0dω = 4π ∫ _∞ 0 r02ρ(r0) 1 r >dr 0 = 4π ∫ _r 0 r02ρ(r0)1 rdr 0_{+ 4π}∫ ∞ r r02ρ(r0)1 r0dr 0 = IntVee(r).

We obtain the matrix of Coulomb potential

Mee=          IntVee(r1) 0 IntVee(r2) . .. 0 IntVee(r399)          399×399 . (2.25)

By (2.19) and (2.20), the integral of IntVee(ri) can be discretized.

2.5.2 Discretization of exchange-correlation potential

The exchange-correlation potential is written as

vxc = εxc[ρ(r)]· ρ(r) δρ(r) = εxc[ρ(r)] + ρ· dεxc[ρ(r)] dρ(r) = εxc[ρ(r)] + ρ· dεxc[ρ(r)] dr · dr dρ(r) = Dif εxc(r), 11

where εxc = εx+ εc is written in section 1.2. Mxc=          Dif εxc(r1) 0 Dif εxc(r2) . .. 0 Dif εxc(r399)          399×399 . (2.26)

By (2.21) and (2.22), the differentiation of Dif εxc(ri) can be discretized.

2.5.3 Discretization of external potential

The external potential is written as vex = −z_r so after the dicretization the

matrix is Mr =          −z r1 0 −z r2 . .. 0 −_rz 399          399×399 , (2.27)

where z is the atomic number.

2.6

Normalization of wave function

The square of a wave function is called the probability density. In quantum mechanics, the probability density can describe the distribution of the elec-trons. In this computation, we have to normalize the wave function to find the wave function that can provide a reasonable probability. The definition of normalization is 1 = ∫ R3 ∫ Pp· drdω,

which the definition of the probability density Pp is

Pp = Ψ∗Ψ = R∗nlRnlYlm∗ Ylm. So that we have 1 = ∫ R3 ∫ Ppdrdω = ∫ _∞ 0 r2R∗_nlRnldr ∫ π o ∫ 2π 0 Y_lm∗ Ylmsinφdθdφ. (2.28) 12

The radial probability density Ppr(r) is defined by

Ppr(r) = r2R∗nlRnl.

Because of the integral [7]

∫ _π

o

∫ _2π

0

Y_lm∗ Ylmsinφdθdφ = 1

and (2.28), we have the result

∫ _∞ 0 Ppr(r)dr = ∫ _∞ 0 r2R∗_nlRnldr = 1. (2.29)

Therefore, we have the conclusion that if we want to normalize a radial wave function Rnl, we just compute the integral

∫ _∞

0

r2R∗_nlRnldr = k, k is a constant,

and divided Rnl by

√

k. Then the new function R0_nl= R_√nl

k is the radial wave

function which has been normalized.

2.7

Self-consistency

Define Mef f be the matrix which signifies the effective potential vef f, and

Mef f = Mee+ Mxc+ Mr. (2.30)

By (2.23), (2.24) and (2.30), the equation (1.18) can be represented as the form

(Mp+ Ml+ Mef f) ~Pnl = enlP~nl,

where the vector ~Pnl is

~ Pnl =          Pnl( ~r1) Pnl( ~r2) .. . Pnl( ~r399)          . It is an eigenvalue problem.

We can determine the correct electronic density and ~Pnl by executing

self-consistency. The table (2.7) can describe the process of self-consistency. In the beginning, build a vector ~ρin to be the initial guess, and we will

get a Hamiltonian Hkh (here Hkh = Mp + Ml + Mef f). Then, computing

the minimal eigenvalue emin of Hkh, and find an eigenvector ~P corresponding

emin. Normalizing ~P and we get the eigenvector ~Pout. Then we can compute

the electronic density ~ρout. Now, if ~ρout = ~ρin it means that we get the right

Pnl, otherwise, generate a new ~ρinas this ~ρout and construct a new Hkh. After

finding the correct electronic density and wave function, we can calculate the ground state energy by using the formula (1.12).

We compute the ground state energy of He (z = 2) to No (z = 102) in the order, and use the preceding atom’s ~ρnlto be the initial guess. For example,

if we want to compute the ground state energy of Mg (z = 12), we will use the ~ρnl of Na (z = 11) to be its initial guess.

out   in



Table (2.7) YES

1. Give an initial guess P_{in nl( )}

2. Compute _in, _{in i}| in nl( ) |2 i occ P f r   



 

3. Construct the Kohn-Sham Hamiltonian H _kh

4. Compute H P_kh_nl _nlP_nl and get P_{out nl( )}

5. Compute _out, _{out i}| out nl( ) |2

i occ P f r   



 

7. Compute the total energy

NO 6. Is_out _in?

16

3 Results

The errors between the computation and the realistic values are less than 15%. The list (3.1) is the comparison of the ideal answer [3] and the answer of this computation with atomic unit, and figure (3.2) is the scale error between the realistic values (Idea data) and the answer of this computat- ion .

The formula of the errors is

errors | i i|, i a b a  

where ai is the data from this computation and bi is the idea data.

Atomic number Idea data

–E(a.u.)

Data in this computation Errors(%)

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

0.500

2.862

7.433

14.57

24.53

37.69

54.40

74.41

99.41

128.5

161.9

199.6

241.9

289

340.7

397.5

459.5

526,8

599.2

676.8

759.7

848.4

942.9

1043

1150

1262

1381

1507

1639

1778

1923

2075

2234

2400

2572

0.500000

2.893138

7.388389

14.67078

25.11337

39.11898

57.24088

79.84583

107.2576

139.8187

173.8986

212.1369

254.2069

300.3957

350.7181

410.3264

464.3733

527.6907

592.0264

659.9070

735.5189

815.8062

900.7615

995.9873

1085.262

1270.632

1384.305

1503.796

1640.438

1761.187

1886.647

2016.884

2151.999

2292.058

2437.159

0.0000 1.0800 0.6000 0.6900 2.3200 3.6500 4.9600 6.8100 7.3200 8.1000 6.9000 5.9100 4.8400 3.7900 2.8600 3.1300 1.0500 0.1700 1.2100 2.5600 3.2900 4.0000 4.6800 4.7200 5.9700 0.6800 0.2400 0.2100 0.0900 0.9500 1.9300 2.8800 3.8100 4.7100 5.5300

18

Atomic number Idea data

-E(a.u.)

Data in this computation Errors(%)

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

2752

2938

3132

3332

3539

3754

3975

4205

4441

4686

4938

5198

5465

5740

6023

6314

6612

6918

7232

7554

7884

8221

8567

8921

9284

9655

10035

10423

10820

11226

11641

12065

12498

12940

13391

2657.920

2809.000

2964.275

3134.383

3310.443

3504.019

3692.796

3887.294

4087.941

4294.852

4523.195

4728.309

4937.590

5155.548

5379.200

5713.260

5953.587

6200.008

6452.556

6697.584

6947.461

7215.920

7535.974

7844.062

8159.394

8482.952

8814.835

9155.146

9476.594

9860.293

10227.65

10602.66

10986.59

11379.53

11781.56

3.5400 4.5900 5.6600 6.3000 6.9000 7.1300 7.6400 8.1700 8.6400 9.1100 9.1700 9.9300 10.6800 11.3400 11.9700 10.5100 11.0600 11.5800 12.0800 12.7900 13.4800 13.9300 13.6800 13.7300 13.7800 13.8200 13.8400 13.8500 14.1800 13.8500 13.8200 13.7900 13.7600 13.7100 13.6600

Atomic number Idea data

-E(a.u.)

Data in this computation Errors(%)

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

13852

14323

14800

15287

15784

16293

16806

17333

17866

18409

18962

19524

20096

20676

21267

21867

22476

23094

23722

24360

25007

25664

26331

27008

27696

28392

29100

29817

30545

31283

32031

32790

12159.12

12542.46

12933.05

13330.92

13736.11

15464.24

15932.53

16442.05

16929.71

17391.06

17868.97

18354.69

18848.17

19349.36

19858.31

20375.08

20871.18

21373.11

21907.81

22450.80

23074.43

23674.40

24320.71

24940.68

25568.91

26166.23

26811.44

27510.10

28178.74

28855.63

29542.39

30239.09

13.9200 14.2000 14.4400 14.6700 14.9100 5.3600 5.4800 5.4200 5.5300 5.8500 6.1200 6.3700 6.6200 6.8600 7.0900 7.3200 7.6900 8.0500 8.2800 8.5000 8.3800 8.4000 8.2700 8.2900 8.3200 8.5100 8.5400 8.3900 8.4000 8.4100 8.4200 8.4400 List (3.1)

References

[1] P.Marder, Micheal. Condensed Matter Physics. A Wiley-Tnterscience Publication, 203-253. 1960.

[2] Gasiorowicz, Stephen. Quantum Physics. John Wiley & Sons, Inc, 197-200. 1995.

[3] Fisher, Forse. Atomic Data. erratum Atom. Data. Nucl. Tables 12 87. 301. 1972.

[4] R. L. Flyrry Jr. Quantum chemistry- An introduction. Englewood Cliffs, NJ: Prentice-Hall, 32-85. 1983.

[5] Parr, Robert G. Density-Function Theory of Atoms and Molecules. New York: Oxford University Press, 1989.

[6] Fiolhais, Carlos, Nogueira. Primer in Density Functional Theory. Berlin; New York: Springer,. 2003.

[7] Menguc, Mustafa, Pinar. Modeling of Radiate Heat transfer in

Multi-didimensional Enclosures using Spherical Harmonics Approximation. A

Bell & Howel information Company, 17. 1985.

[8] Gilbarg David. Elliptic partial differential equations of second order. New York: Springer, 2001.

[9] Richard, Courant and David, Hilbert. Methods of mathematical Physics. 1962.

[10] Dahl, Jens Peder ed. Local Density Approximation in Quantum

Chem-istry and Solid State Physics. New York: Plenum, 1984.

Appendix

The source file of this computation module global

implicit none

integer ,parameter ::d=400 integer ,parameter ::sz=399 integer i,j,Bool,k,times real ,parameter ::pi=3.1416 real ,parameter ::bdd =10.0 real Yini(sz),Nini(sz),Pin(sz),Rin(sz),Yin(sz),ABNd(sz) real Nin(sz),Pout(sz),Rout(sz),Yout(sz),Nout(sz),R(sz),N(sz),Nd(sz) real Sch(sz,sz),L1(sz,sz),L2(sz,sz) real Veff(sz,sz),Exc(sz,sz),Eee(sz,sz),Rv(sz,sz) real s,Ncin,Ncout real A_2nd(sz),B_2nd(sz),A_3rd(sz),B_3rd(sz) real A_4th(sz),B_4th(sz),A_5th(sz),B_5th(sz) real A_6th(sz),B_6th(sz),A_7th(sz),B_7th(sz) integer k_2nd,k_3rd,k_4th,k_5th,k_6th,k_7th real Pin_2s(sz),Pout_2s(sz),H_2s(sz),Rout_2s(sz) real Pout_3s(sz),Pin_3s(sz),Rout_3s(sz),H_3s(sz) real Pin_2p(sz),Pout_2p(sz),H_2p(sz),Rout_2p(sz) real H_3p(sz),Pout_3p(sz),Pin_3p(sz),Rout_3p(sz) real Pin_4s(sz),Pout_4s(sz),H_4s(sz),Rout_4s(sz) real Pin_3d(sz),Pout_3d(sz),H_3d(sz),Rout_3d(sz) real Pin_4p(sz),Pout_4p(sz),H_4p(sz),Rout_4p(sz) real Pin_5s(sz),Pout_5s(sz),H_5s(sz),Rout_5s(sz) real Pin_4d(sz),Pout_4d(sz),H_4d(sz),Rout_4d(sz) real Pin_5p(sz),Pout_5p(sz),H_5p(sz),Rout_5p(sz) real Pin_6s(sz),Pout_6s(sz),H_6s(sz),Rout_6s(sz) real Pin_5d(sz),Pout_5d(sz),H_5d(sz),Rout_5d(sz) real Pin_4f(sz),Pout_4f(sz),H_4f(sz) real Rout_4f(sz) real Pin_6p(sz),Pout_6p(sz),H_6p(sz),Rout_6p(sz) real Pin_7s(sz),Pout_7s(sz),H_7s(sz),Rout_7s(sz) real Pin_6d(sz),Pout_6d(sz),H_6d(sz),Rout_6d(sz)

real Pin_5f(sz),Pout_5f(sz),H_5f(sz),Rout_5f(sz) real Pin_7p(sz),Pout_7p(sz),H_7p(sz),Rout_7p(sz) !sch 變數 real H(sz),eigenvalue(sz),eigenvector(sz,sz) !L1 變數 real hg !L2 變數 !Veff 變數 !Ec 變數 real Ec(sz,sz)

real ,parameter ::Ge=0.1423 real ,parameter ::Be1=1.0529 real ,parameter ::Be2=0.3334 real Ech(sz,sz)

real ,parameter ::Aec=0.0311 real ,parameter ::Bec=-0.048 real ,parameter ::Cec=0.002 real ,parameter ::Dec=-0.0116 !Ex 變數 real Rs(sz),Ex(sz,sz) !Ediff 變數 real Addxc(sz,sz),Edr(sz,sz),Ndr(sz,sz),Ediff(sz,sz) !Exc 變數 !intg1 變數

!real ,external ::Intg1 !intg1 變數

!real ,external ::Intg2 !Eee 變數 !Exc 變數 !Energy 的變數 real total real lPout_1sl(sz),lPout_2sl(sz),lPout_3sl(sz),lPout_4sl(sz),lPout_5sl(sz) ,lPout_6sl(sz),lPout_7sl(sz) real lPout_2pl(sz),lPout_3pl(sz),lPout_4pl(sz),lPout_5pl(sz),lPout_6pl(sz) real lPout_3dl(sz),lPout_4dl(sz),lPout_5dl(sz),lPout_6dl(sz) 22

real lPout_4fl(sz),lPout_5fl(sz) real Hn(sz),Hn_2s(sz),Hn_3s(sz),Hn_4s(sz),Hn_5s(sz),Hn_6s(sz),Hn_7s(sz) real Hn_2p(sz),Hn_3p(sz),Hn_4p(sz),Hn_5p(sz),Hn_6p(sz) real Hn_3d(sz),Hn_4d(sz),Hn_5d(sz),Hn_6d(sz) real Hn_4f(sz),Hn_5f(sz) end module !--- Program main use IMSL use global implicit none call self_he() open(unit=11,file='Reslet1.txt') open(unit=12,file='eigenvalueofHe.txt') call total_energy(2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"He(2)" write(*,*)times write(*,*)total write(11,*)"He(2)" write(11,*)times write(11,*)total !do i=1,sz !write(12,*)eigenvalue(i) !end do call Li3() call Be4() call B5() call C6() call N7() call O8() call F9() call Ne10()

call Na11() call Mg12() call Al13() call Si14() call P15() call S16() call Cl17() call Ar18() call K19() call Ca20() call Sc21() call Ti22() call V23() call Cr24() call Mn25() call Fe26() call Co27() call Ni28() call Cu29() call Zn30() call Ga31() call Ge32() call As33() call Se34() call Br35() call Kr36() call Rb37() call Sr38() call Y39() call Zr40() call Nb41() call Mo42() call Tc43() call Ru44() call Rh45() call Pb46() call Ag47() call Cd48() 24

call In49() call Sn50() call Sb51() call Te52() call I53() call Xe54() call Cs55() call Ba56() call La57() call Ce58() call Pr59() call Nd60() call Pm61() call Sm62() call Eu63() call Gd64() call Tb65() call Dy66() call Ho67() call Er68() call Tm69() call Yb70() call Lu71() call Hf72() call Ta73() call W74() call Re75() call Os76() call Ir77() call Pt78() call Au79() call Hg80() call Tl81() call Pb82() call Bi83() call Po84() call At85() call Rn86()

call Fr87() call Ra88() call Ac89() call Th90() call Pa91() call U92() call Np93() call Pu94() call Am95() call Cm96() call Bk97() call Cf98() call Es99() call Fm100() call Md101() call No102() call Lr103() call Rf104() call Db105() call Sg106() call Bh107() end !--- subroutine subSch_7s() !矩陣合 use global use IMSL implicit none

real ,external :: norl do i=1,sz do j=1,sz Sch(i,j)=L1(i,j)+L2(i,j)+Veff(i,j) end do end do eigenvalue=eig(Sch,V=eigenvector) !為了存取特徵值的最小值位置 open(unit=10,file="7s_1st.txt") 26

write(10,*)minloc(eigenvalue) rewind(10) read(10,*)k !--- do i=1,sz H(i)=eigenvector(i,k) end do do i=1,sz Hn(i)=H(i)/norl(R,H) end do do i=1,sz A_2nd(i)=0 end do A_2nd(k)=100000 do i=1,sz B_2nd(i)=A_2nd(i)+eigenvalue(i) end do open(unit=10,file="7s_2nd.txt") write(10,*)minloc(B_2nd) rewind(10) read(10,*)k_2nd do i=1,sz H_2s(i)=eigenvector(i,k_2nd) end do do i=1,sz Hn_2s(i)=H_2s(i)/norl(R,H_2s) end do !--- do i=1,sz A_3rd(i)=0 end do A_3rd(k_2nd)=100000 do i=1,sz B_3rd(i)=B_2nd(i)+A_3rd(i) end do open(unit=10,file="7s_3rd.txt") write(10,*)minloc(B_3rd) rewind(10)

read(10,*)k_3rd do i=1,sz H_3s(i)=eigenvector(i,k_3rd) end do do i=1,sz Hn_3s(i)=H_3s(i)/norl(R,H_3s) end do do i=1,sz A_4th(i)=0 end do A_4th(k_3rd)=100000 do i=1,sz B_4th(i)=B_3rd(i)+A_4th(i) end do open(unit=10,file="7s_4th.txt") write(10,*)minloc(B_4th) rewind(10) read(10,*)k_4th do i=1,sz H_4s(i)=eigenvector(i,k_4th) end do do i=1,sz Hn_4s(i)=H_4s(i)/norl(R,H_4s) end do do i=1,sz A_5th(i)=0 end do A_5th(k_4th)=100000 do i=1,sz B_5th(i)=B_4th(i)+A_5th(i) end do open(unit=10,file="7s_5th.txt") write(10,*)minloc(B_5th) rewind(10) read(10,*)k_5th do i=1,sz H_5s(i)=eigenvector(i,k_5th) end do 28

do i=1,sz Hn_5s(i)=H_5s(i)/norl(R,H_5s) end do do i=1,sz A_6th(i)=0 end do A_6th(k_5th)=100000 do i=1,sz B_6th(i)=B_5th(i)+A_6th(i) end do open(unit=10,file="7s_6th.txt") write(10,*)minloc(B_6th) rewind(10) read(10,*)k_6th do i=1,sz H_6s(i)=eigenvector(i,k_6th) end do do i=1,sz Hn_6s(i)=H_6s(i)/norl(R,H_6s) end do do i=1,sz A_7th(i)=0 end do A_7th(k_6th)=100000 do i=1,sz B_7th(i)=B_6th(i)+A_7th(i) end do open(unit=10,file="7s_7th.txt") write(10,*)minloc(B_7th) rewind(10) read(10,*)k_7th do i=1,sz H_7s(i)=eigenvector(i,k_7th) end do do i=1,sz Hn_7s(i)=H_7s(i)/norl(R,H_7s) end do

end

subroutine subSch_6p() !矩陣合 use global

use IMSL implicit none

real ,external :: norl do i=1,sz do j=1,sz Sch(i,j)=L1(i,j)+L2(i,j)+Veff(i,j) end do end do eigenvalue=eig(Sch,V=eigenvector) !為了存取特徵值的最小值位置 open(unit=10,file="6p_1st.txt") write(10,*)minloc(eigenvalue) rewind(10) read(10,*)k !--- do i=1,sz H_2p(i)=eigenvector(i,k) end do do i=1,sz Hn_2p(i)=H_2p(i)/norl(R,H_2p) end do do i=1,sz A_2nd(i)=0 end do A_2nd(k)=100000 do i=1,sz B_2nd(i)=A_2nd(i)+eigenvalue(i) end do open(unit=10,file="6P_2nd.txt") write(10,*)minloc(B_2nd) rewind(10) read(10,*)k_2nd do i=1,sz H_3p(i)=eigenvector(i,k_2nd) end do 30

do i=1,sz Hn_3p(i)=H_3p(i)/norl(R,H_3p) end do do i=1,sz A_3rd(i)=0 end do A_3rd(k_2nd)=100000 do i=1,sz B_3rd(i)=A_3rd(i)+B_2nd(i) end do open(unit=10,file="6P_3rd.txt") write(10,*)minloc(B_3rd) rewind(10) read(10,*)k_3rd do i=1,sz H_4p(i)=eigenvector(i,k_3rd) end do do i=1,sz Hn_4p(i)=H_4p(i)/norl(R,H_4p) end do !--- do i=1,sz A_4th(i)=0 end do A_4th(k_3rd)=100000 do i=1,sz B_4th(i)=A_4th(i)+B_3rd(i) end do open(unit=10,file="6P_4th.txt") write(10,*)minloc(B_4th) rewind(10) read(10,*)k_4th do i=1,sz H_5p(i)=eigenvector(i,k_4th) end do do i=1,sz Hn_5p(i)=H_5p(i)/norl(R,H_5p) end do

do i=1,sz A_5th(i)=0 end do A_5th(k_4th)=100000 do i=1,sz B_5th(i)=A_5th(i)+B_4th(i) end do open(unit=10,file="6P_5th.txt") write(10,*)minloc(B_5th) rewind(10) read(10,*)k_5th do i=1,sz H_6p(i)=eigenvector(i,k_5th) end do do i=1,sz Hn_6p(i)=H_6p(i)/norl(R,H_6p) end do end subroutine subSch_6d() !矩陣合 use global use IMSL implicit none

real ,external ::norl do i=1,sz do j=1,sz Sch(i,j)=L1(i,j)+L2(i,j)+Veff(i,j) end do end do eigenvalue=eig(Sch,V=eigenvector) !為了存取特徵值的最小值位置 open(unit=10,file="6d_1st.txt") write(10,*)minloc(eigenvalue) rewind(10) read(10,*)k !--- do i=1,sz H_3d(i)=eigenvector(i,k) end do 32

do i=1,sz Hn_3d(i)=H_3d(i)/norl(R,H_3d) end do do i=1,sz A_2nd(i)=0 end do A_2nd(k)=100000 do i=1,sz B_2nd(i)=A_2nd(i)+eigenvalue(i) end do open(unit=10,file="6d_2nd.txt") write(10,*)minloc(B_2nd) rewind(10) read(10,*)k_2nd do i=1,sz H_4d(i)=eigenvector(i,k_2nd) end do do i=1,sz Hn_4d(i)=H_4d(i)/norl(R,H_4d) end do do i=1,sz A_3rd(i)=0 end do A_3rd(k_2nd)=100000 do i=1,sz B_3rd(i)=A_3rd(i)+B_2nd(i) end do open(unit=10,file="6d_3rd.txt") write(10,*)minloc(B_3rd) rewind(10) read(10,*)k_3rd do i=1,sz H_5d(i)=eigenvector(i,k_3rd) end do do i=1,sz Hn_5d(i)=H_5d(i)/norl(R,H_5d) end do do i=1,sz

A_4th(i)=0 end do A_4th(k_3rd)=100000 do i=1,sz B_4th(i)=A_4th(i)+B_3rd(i) end do open(unit=10,file="6d_4th.txt") write(10,*)minloc(B_4th) rewind(10) read(10,*)k_4th do i=1,sz H_6d(i)=eigenvector(i,k_4th) end do do i=1,sz Hn_6d(i)=H_6d(i)/norl(R,H_6d) end do end subroutine subSch_5f() !矩陣合 use global use IMSL implicit none

real ,external :: norl do i=1,sz do j=1,sz Sch(i,j)=L1(i,j)+L2(i,j)+Veff(i,j) end do end do eigenvalue=eig(Sch,V=eigenvector) !為了存取特徵值的最小值位置 open(unit=10,file="5f_1st.txt") write(10,*)minloc(eigenvalue) rewind(10) read(10,*)k !--- do i=1,sz H_4f(i)=eigenvector(i,k) end do 34

do i=1,sz Hn_4f(i)=H_4f(i)/norl(R,H_4f) end do do i=1,sz A_2nd(i)=0 end do A_2nd(k)=100000 do i=1,sz B_2nd(i)=A_2nd(i)+eigenvalue(i) end do open(unit=10,file="5f_2nd.txt") write(10,*)minloc(B_2nd) rewind(10) read(10,*)k_2nd do i=1,sz H_5f(i)=eigenvector(i,k_2nd) end do do i=1,sz Hn_5f(i)=H_5f(i)/norl(R,H_5f) end do end subroutine subL1() !兩次倒數用泰勒展開之後逼近的結果 use IMSL use global implicit none L1=eye(sz) hg=bdd/d do i=1,sz L1(i,i)=1.0/(hg**2.0) end do do i=2,sz L1(i,i-1)=-1.0/(2.0*(hg**2.0)) end do do i=1,sz-1 L1(i,i+1)=-1.0/(2.0*(hg**2.0)) end do end

subroutine subL2(l) !l(l+1)/2r^2 那項 use IMSL use global implicit none integer l L2=eye(sz) hg=bdd/d do i=1,sz L2(i,i)=(l*(l+1))/(2.0*((hg*i)**2.0)) end do end

subroutine subVeff(z) !Veff=Vee+Vxc+V 那項(總合) use IMSL use global implicit none integer z hg=bdd/d !-z/r(Vr 那項) Rv=eye(sz) do i=1,sz Rv(i,i)=-z/(hg*i) end do do i=1,sz do j=1,sz Veff(i,j)=Eee(i,j)+Exc(i,j)+Rv(i,j) end do end do end

subroutine subExc() !Vex=Exc+n*dExc/dn(總合) use global use IMSL implicit none do j=1,sz do i=1,sz If (Rs(i)>=1) then Exc(i,j)=Ex(i,j)+Ec(i,j)+Ediff(i,j) else Exc(i,j)=Ex(i,j)+Ech(i,j)+Ediff(i,j) 36

end if end do end do end subroutine subEx() use IMSL use global implicit none Ex=eye(sz) do i=1,sz rs(i)=(3.0/(4.0*pi*N(i)))**(1.0/3.0) end do do i=1,sz Ex(i,i)=(-0.458)/(rs(i)) end do end subroutine subEc() use IMSL use global implicit none Ec=eye(sz) do i=1,sz Ec(i,i)=-Ge/(1.0+Be1*(Rs(i)**(0.5))+Be2*Rs(i)) end do end subroutine subEch() use IMSL use global implicit none Ech=eye(sz) do i=1,sz Ech(i,i)=Aec*LOG(Rs(i))+Bec+Cec*Rs(i)*LOG(Rs(i))+Dec*Rs(i) end do end subroutine subEdiff() !用數值微分 use global use IMSL implicit none

Ediff=eye(sz) do i=1,sz R(i)=(bdd/d)*i end do Edr=eye(sz) Ndr=eye(sz) do j=1,sz do i=1,sz if (Rs(i)>=1) then Addxc(i,j)=Ex(i,j)+Ec(i,j) !其實就是 Exc else Addxc(i,j)=Ex(i,j)+Ech(i,j) end if end do end do Edr(1,1)=(Addxc(1,1)-Addxc(2,2))/(R(1)-R(2)) Edr(sz,sz)=(Addxc(sz-1,sz-1)-Addxc(sz,sz))/(R(sz-1)-R(sz)) do i=2,sz-1 Edr(i,i)=(Addxc(i+1,i+1)-Addxc(i-1,i-1))/(R(i+1)-R(i-1)) end do Ndr(1,1)=(N(1)-N(2))/(R(1)-R(2)) Ndr(sz,sz)=(N(sz-1)-N(sz))/(R(sz-1)-R(sz)) do i=2,sz-1 Ndr(i,i)=(N(i+1)-N(i-1))/(R(i+1)-R(i-1)) end do do i=1,sz Ediff(i,i)=Edr(i,i)*(1/Ndr(i,i))*N(i) end do end subroutine subEee() use Global use IMSL implicit none real ,external::Intg1 real ,external::Intg2 Eee=eye(sz) do i=1,sz R(i)=(bdd/d)*i 38

end do do i=1,sz

Eee(i,i)=Intg1(i,n)+Intg2(i,n) end do

end

real function Intg1(a,n) !Eee 的數值積分 use IMSL implicit none integer ,parameter::d=400 integer ,parameter::sz=399 real ,parameter::bdd=10.0 real N(sz) real hg,Int

real ,parameter ::pi=3.1416 integer i,a Int=0.0 hg=bdd/d do i=1,a-1 Int=Int+((((hg*i)**2.0)*N(i)/(a*hg))+(((hg*(i+1))**2.0)*N(i+1)/(a*hg) )) *hg*0.5*4*pi end do Intg1=Int end function

real function Intg2(b,n) use IMSL

implicit none

integer ,parameter ::d=400 integer ,parameter ::sz=399 real ,parameter::bdd=10.0 real ,parameter ::pi=3.1416 real N(sz) real hg,Int integer b,i Int=0.0 hg=bdd/d do i=b,sz-1 Int=Int+((((hg*i)**2.0)*N(i))/(hg*i)+(((hg*(i+1))**2.0)*N(i+1))/(hg*( i+1))) *hg*0.5*4*pi

end do

Intg2=Int !+( (((hg*sz)**2.0)*N(sz))/(hg*sz) )*0.5*hg end function

subroutine subProcess_6p(l,z)!所有的 subroutine 都呼叫一次 use IMSL use global implicit none integer l,z call subL1() call subL2(l) call subEee() call subEx() call subEc() call subEdiff() call subExc() call subVeff(z) call subSch_6p() end

subroutine subProcess_7s(l,z)!所有的 subroutine 都呼叫一次 use IMSL use global implicit none integer l,z call subL1() call subL2(l) call subEee() call subEx() call subEc() call subEdiff() call subExc() call subVeff(z) call subSch_7s() end

subroutine subProcess_6d(l,z)!所有的 subroutine 都呼叫一次 use IMSL

use global implicit none integer l,z

call subL1() call subL2(l) call subEee() call subEx() call subEc() call subEdiff() call subExc() call subVeff(z) call subSch_6d() end

subroutine subProcess_5f(l,z)!所有的 subroutine 都呼叫一次 use IMSL use global implicit none integer l,z call subL1() call subL2(l) call subEee() call subEx() call subEc() call subEdiff() call subExc() call subVeff(z) call subSch_5f() end !-Total Energy--- subroutine Total_Energy(z,N1s,N2s,N2p,N3s,N3p,N4s,N3d,N4p,N5s,N4d,N5p,N6s,N5d,N4 f,N6p,N7s,N6d,N5f)!4d 的 energy use IMSL use global implicit none integer z,N1s,N2s,N3s,N4s,N5s,N6s,N7s,N2p,N3p,N4p,N5p,N6p,N3d,N4d,N5d,N6d,N4f ,N5f

real ,external :: Env real ,external :: Enee real ,external :: Enxc

real ,external :: Intk1 real ,external :: Intk2 real ,external :: norl real Enkk do i=1,sz Rout(i)=Pout(i)/R(i) Rout_2s(i)=Pout_2s(i)/R(i) Rout_3s(i)=Pout_3s(i)/R(i) Rout_3p(i)=Pout_3p(i)/R(i) Rout_4s(i)=Pout_4s(i)/R(i) Rout_3d(i)=Pout_3d(i)/R(i) Rout_4p(i)=Pout_4p(i)/R(i) Rout_5s(i)=Pout_5s(i)/R(i) Rout_4d(i)=Pout_4d(i)/R(i) Rout_5p(i)=Pout_5p(i)/R(i) Rout_6s(i)=Pout_6s(i)/R(i) Rout_4f(i)=Pout_4f(i)/R(i) Rout_5d(i)=Pout_5d(i)/R(i) Rout_6p(i)=Pout_6p(i)/R(i) Rout_7s(i)=Pout_7s(i)/R(i) Rout_5f(i)=Pout_5f(i)/R(i) Rout_6d(i)=Pout_6d(i)/R(i) end do Enkk = N1s*( Intk1(Rout,R)+Intk2(Rout,R,0) )+N2s*( Intk1(Rout_2s,R)+Intk2(Ro ut_2s,R,0) )+N2p*( Intk1(Rout_2p,R)+Intk2(Rout_2p,R,1) )+N3s*( Intk1( Rout_3s,R)+Intk2(Rout_3s,R,0) )+N3p*( Intk1(Rout_3p,R)+Intk2(Rout_3p, R,1) )+N4s*( Intk1(Rout_4s,R)+Intk2(Rout_4s,R,0) )+N3d*( Intk1(Rout_3 d,R)+Intk2(Rout_3d,R,2) )+N4p*( Intk1(Rout_4p,R)+Intk2(Rout_4p,R,1) ) +N5s*( Intk1(Rout_5s,R)+Intk2(Rout_5s,R,0) )+N4d*( Intk1(Rout_4d,R)+I ntk2(Rout_4d,R,2) )+N5p*( Intk1(Rout_5p,R)+Intk2(Rout_5p,R,1) )+N6s*( Intk1(Rout_6s,R)+Intk2(Rout_6s,R,0))+N5d*(Intk1(Rout_5d,R)+Intk2(Rout _5d,R,2))+N4f*(Intk1(Rout_4f,R)+Intk2(Rout_4f,R,3))+N6p*(Intk1(Rout_6 p,R)+Intk2(Rout_6p,R,1))+N7s*( Intk1(Rout_7s,R)+Intk2(Rout_7s,R,0) )+ N6d*( Intk1(Rout_6d,R)+Intk2(Rout_6d,R,2) )+N5f*( Intk1(Rout_5f,R)+In tk2(Rout_5f,R,3) ) total = Enkk+Env(Nout,R,z)+Enxc(Nout,R,Addxc)+Enee(Nout,R,z) end 42

real function Env(n,r,z)!(外位能) implicit none

integer i,z

real ,parameter ::pi=3.1416 integer ,parameter ::d=400 integer ,parameter ::sz=399 real bdd real n(sz),r(sz) Env=0.0 bdd=10.0 do i=1,sz-1

Env = Env-( z*r(i)*n(i)+z*r(i+1)*n(i+1) )*4*pi*0.5*(bdd/d) end do

end function

real function Enxc(n,r,x)!Exc implicit none

integer i

real ,parameter ::pi=3.1416 integer ,parameter ::d=400 integer ,parameter ::sz=399 real hg,bdd real n(sz),r(sz),x(sz,sz) bdd=10.0 hg=bdd/d Enxc=0.0 do i=1,sz-1 Enxc=Enxc+4*pi*( x(i,i)*n(i)*(r(i)**2)+x(i+1,i+1)*n(i+1)*(r(i+1)**2) )*0.5*(bdd/d) end do end function

real function Enee(n,r,z)!Eee implicit none

integer i,j,z

real ,parameter ::pi=3.1416 integer ,parameter ::d=400 integer ,parameter ::sz=399

real n(sz),r(sz),bdd bdd=10.0

Enee=0.0 do i=1,sz do j=1,sz

Enee=Enee+0.5*( (( r(i)*r(j) )**2)/( max(i,j)*(bdd/d) ) )*n(i)*n(j)*( bdd/d)*(bdd/d)*4*pi*4*pi

end do end do end

real function Intk1(Rc,r) implicit none

integer i,j

real ,parameter ::pi=3.1416 integer ,parameter ::d=400 integer ,parameter ::sz=399 real hg real n(sz),r(sz),bdd,Rc(sz),inti1 bdd=10.0 inti1=0.0 hg=bdd/d do i=2,sz-2

Inti1=Inti1+(( -0.5*Rc(i)*r(i)*( (Rc(i+1)-Rc(i-1))/(r(i+1)-r(i-1))+(R c(i+1)-2*Rc(i)+Rc(i-1))/(hg**2) ) )+( -0.5*Rc(i+1)*r(i+1)*( (Rc(i+2) -Rc(i))/(r(i+2)-r(i))+(Rc(i+2)-2*Rc(i+1)+Rc(i))/(hg**2) ) ))*hg*0.5 !( Rc(i)*( (-1.0/r(i))*( (Rc(i-1)-Rc(i+1))/(R(i-1)-R(i+1)) ) )-0.5 *( (Rc(i+1)-2*Rc(i)+Rc(i-1))/((bdd/d)**2) )*hg*0.5 ) +

( Rc(i+1)*( (-1.0/r(i+1))*( (Rc(i)-Rc(i+2))/(R(i)-R(i+2)) ) )-0.5*( (Rc(i+2)-2*Rc(i+1)+Rc(i))/((bdd/d)**2))*hg*0.5 )

end do Intk1=inti1 end function

real function Intk2(Rc,r,l) implicit none

integer i,j,l

real ,parameter ::pi=3.1416 integer ,parameter ::d=400

integer ,parameter ::sz=399 real hg,Intk real n(sz),r(sz),bdd,Rc(sz) bdd=10.0 intk=0.0 hg=bdd/d do i=1,sz-1

Intk=Intk+( ( (Rc(i)**2)*l*(l+1)*0.5 )+( (Rc(i+1)**2)*l*(l+1)*0.5 ) )*0.5*hg end do Intk2=Intk end function subroutine self_he() use IMSL use global implicit none

real ,external :: norl ! bool=1 !自洽迴圈用的值 hg=bdd/d !空隙 do i=1,sz R(i)=hg*i !R 向量 domain 的值 end do !以下為初始值 do i=1,sz Yini(i)=exp(-R(i)/0.53)!*2.0*((2.0/0.53)**(1.5)) end do do i=1,sz Nini(i)=2.0*((Yini(i))**2.0) end do !---初始值帶入矩陣 do i=1,sz N(i)=Nini(i) end do call subProcess_7s(0,2) do i=1,sz Pout(i)=Hn(i) end do

do i=1,sz Nout(i)=2.0*((Pout(i)/R(i))**2.0) end do !(以上)----初始直的結果-- !-給 s 初值-- s=10 times=0 !---自洽---He do while (s>=10**(-3.0)) !誤差值在 10**-6 如果寫 s=0 跑不完 times=times+1 do i=1,sz N(i)=Nout(i) end do call subProcess_7s(0,2) do i=1,sz Pin(i)=Hn(i) end do do i=1,sz Nin(i)=2.0*((Pin(i)/R(i))**2.0) end do do i=1,sz Nd(i)=Nin(i)-Nout(i) end do do i=1,sz ABNd(i)=ABS(Nd(i)) end do s=maxval(ABNd) !write(*,*)times !write(*,*)s !--- do i=1,sz Nout(i)=Nin(i) end do end do end subroutine Li3() use IMSL 46

use global implicit none !---自洽---Li call sub_self(3,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(3,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Li(3)" write(*,*)times write(*,*)total write(11,*)"Li(3)" write(11,*)times write(11,*)total end subroutine Be4() use IMSL use global implicit none !---自洽---Be call sub_self(4,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(4,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Be(4)" write(*,*)times write(*,*)total write(11,*)"Be(4)" write(11,*)times write(11,*)total end subroutine B5() use IMSL use global implicit none !----B 自洽 --- call sub_self(5,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(5,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"B(5)" write(*,*)times

write(*,*)total write(11,*)"B(5)" write(11,*)times write(11,*)total end subroutine C6() use IMSL use global implicit none !---C call sub_self(6,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(6,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"C(6)" write(*,*)times write(*,*)total write(11,*)"C(6)" write(11,*)times write(11,*)total end !---N---- subroutine N7() use IMSL use global implicit none call sub_self(7,2,2,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(7,2,2,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"N(7)" write(*,*)times write(*,*)total write(11,*)"N(7)" write(11,*)times write(11,*)total end subroutine O8() 48

use IMSL use global implicit none !---O---- call sub_self(8,2,2,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(8,2,2,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"O(8)" write(*,*)times write(*,*)total write(11,*)"O(8)" write(11,*)times write(11,*)total end subroutine F9() use IMSL use global implicit none !----F--- call sub_self(9,2,2,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(9,2,2,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"F(9)" write(*,*)times write(*,*)total write(11,*)"F(9)" write(11,*)times write(11,*)total end subroutine Ne10() use IMSL use global implicit none !---Ne--- call sub_self(10,2,2,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(10,2,2,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Ne(10)"

write(*,*)times write(*,*)total write(11,*)"Ne(10)" write(11,*)times write(11,*)total end subroutine Na11() use IMSL use global implicit none !---Na--- call sub_self(11,2,2,6,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(11,2,2,6,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Na(11)" write(*,*)times write(*,*)total write(11,*)"Na(11)" write(11,*)times write(11,*)total end subroutine Mg12() use IMSL use global implicit none !---Mg--- call sub_self(12,2,2,6,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(12,2,2,6,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Mg(12)" write(*,*)times write(*,*)total write(11,*)"Mg(12)" write(11,*)times write(11,*)total end subroutine AL13() 50

use IMSL use global implicit none !----AL---- call sub_self(13,2,2,6,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(13,2,2,6,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"AL(13)" write(*,*)times write(*,*)total write(11,*)"AL(13)" write(11,*)times write(11,*)total end subroutine SI14() use IMSL use global implicit none !----SI---- call sub_self(14,2,2,6,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(14,2,2,6,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"SI(14)" write(*,*)times write(*,*)total write(11,*)"SI(14)" write(11,*)times write(11,*)total end subroutine P15() use IMSL use global implicit none !----P---- call sub_self(15,2,2,6,2,3,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(15,2,2,6,2,3,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"P(15)" write(*,*)times

write(*,*)total write(11,*)"P(15)" write(11,*)times write(11,*)total end subroutine S16() use IMSL use global implicit none !----S---- call sub_self(16,2,2,6,2,4,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(16,2,6,2,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"S(16)" write(*,*)times write(*,*)total write(11,*)"S(16)" write(11,*)times write(11,*)total end subroutine CL17() use IMSL use global implicit none !----CL---- call sub_self(17,2,2,6,2,5,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(17,2,2,6,2,5,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"CL(17)" write(*,*)times write(*,*)total write(11,*)"CL(17)" write(11,*)times write(11,*)total end subroutine Ar18() use IMSL 52

use global implicit none !----Ar---- call sub_self(18,2,2,6,2,6,0,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(18,2,2,6,2,6,0,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Ar(18)" write(*,*)times write(*,*)total write(11,*)"Ar(18)" write(11,*)times write(11,*)total end subroutine K19() use IMSL use global implicit none !---K---- call sub_self(19,2,2,6,2,6,1,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(19,2,2,6,2,6,1,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"K(19)" write(*,*)times write(*,*)total write(11,*)"K(19)" write(11,*)times write(11,*)total end subroutine Ca20() use IMSL use global implicit none !---Ca--- call sub_self(20,2,2,6,2,6,2,0,0,0,0,0,0,0,0,0,0,0,0) call total_energy(20,2,2,6,2,6,2,0,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Ca(20)" write(*,*)times write(*,*)total

write(11,*)"Ca(20)" write(11,*)times write(11,*)total end subroutine Sc21() use IMSL use global implicit none !---Sc--- call sub_self(21,2,2,6,2,6,2,1,0,0,0,0,0,0,0,0,0,0,0) call total_energy(21,2,2,6,2,6,2,1,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Sc(21)" write(*,*)times write(*,*)total write(11,*)"Sc(21)" write(11,*)times write(11,*)total end subroutine Ti22() use IMSL use global implicit none !---Ti--- call sub_self(22,2,2,6,2,6,2,2,0,0,0,0,0,0,0,0,0,0,0) call total_energy(22,2,2,6,2,6,2,2,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Ti(22)" write(*,*)times write(*,*)total write(11,*)"Ti(22)" write(11,*)times write(11,*)total end subroutine V23() use IMSL use global 54

implicit none !---V--- call sub_self(23,2,2,6,2,6,2,3,0,0,0,0,0,0,0,0,0,0,0) call total_energy(23,2,2,6,2,6,2,3,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"V(23)" write(*,*)times write(*,*)total write(11,*)"V(23)" write(11,*)times write(11,*)total end subroutine Cr24() use IMSL use global implicit none !---Cr---- call sub_self(24,2,2,6,2,6,1,5,0,0,0,0,0,0,0,0,0,0,0) call total_energy(24,2,2,6,2,6,1,5,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Cr(24)" write(*,*)times write(*,*)total write(11,*)"Cr(24)" write(11,*)times write(11,*)total end subroutine Mn25() use IMSL use global implicit none !---Mn--- call sub_self(25,2,2,6,2,6,2,5,0,0,0,0,0,0,0,0,0,0,0) call total_energy(25,2,2,6,2,6,2,5,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Mn(25)" write(*,*)times write(*,*)total write(11,*)"Mn(25)"

write(11,*)times write(11,*)total end subroutine Fe26() use IMSL use global implicit none !---Fe--- call sub_self(26,2,2,6,2,6,2,6,0,0,0,0,0,0,0,0,0,0,0) call total_energy(26,2,2,6,2,6,2,6,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Fe(26)" write(*,*)times write(*,*)total write(11,*)"Fe(26)" write(11,*)times write(11,*)total end subroutine Co27() use IMSL use global implicit none !---Co--- call sub_self(27,2,2,6,2,6,2,7,0,0,0,0,0,0,0,0,0,0,0) call total_energy(27,2,2,6,2,6,2,7,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Co(27)" write(*,*)times write(*,*)total write(11,*)"Co(27)" write(11,*)times write(11,*)total end subroutine Ni28() use IMSL use global implicit none 56

!---Ni--- call sub_self(28,2,2,6,2,6,2,8,0,0,0,0,0,0,0,0,0,0,0) call total_energy(28,2,2,6,2,6,2,8,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Ni(28)" write(*,*)times write(*,*)total write(11,*)"Ni(28)" write(11,*)times write(11,*)total end subroutine Cu29() use IMSL use global implicit none !---Cu--- call sub_self(29,2,2,6,2,6,1,10,0,0,0,0,0,0,0,0,0,0,0) call total_energy(29,2,2,6,2,6,1,10,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Cu(29)" write(*,*)times write(*,*)total write(11,*)"Cu(29)" write(11,*)times write(11,*)total end subroutine Zn30() use IMSL use global implicit none !---Zn--- call sub_self(30,2,2,6,2,6,2,10,0,0,0,0,0,0,0,0,0,0,0) call total_energy(30,2,2,6,2,6,2,10,0,0,0,0,0,0,0,0,0,0,0) write(*,*)"Zn(30)" write(*,*)times write(*,*)total write(11,*)"Zn(30)" write(11,*)times

write(11,*)total end subroutine Ga31() use IMSL use global implicit none !---Ga--- call sub_self(31,2,2,6,2,6,2,10,1,0,0,0,0,0,0,0,0,0,0) call total_energy(31,2,2,6,2,6,2,10,1,0,0,0,0,0,0,0,0,0,0) write(*,*)"Ga(31)" write(*,*)times write(*,*)total write(11,*)"Ga(31)" write(11,*)times write(11,*)total end subroutine Ge32() use IMSL use global implicit none !---Ge--- call sub_self(32,2,2,6,2,6,2,10,2,0,0,0,0,0,0,0,0,0,0) call total_energy(32,2,2,6,2,6,2,10,2,0,0,0,0,0,0,0,0,0,0) write(*,*)"Ge(32)" write(*,*)times write(*,*)total write(11,*)"Ge(32)" write(11,*)times write(11,*)total end subroutine As33() use IMSL use global implicit none !---As--- 58

call sub_self(33,2,2,6,2,6,2,10,3,0,0,0,0,0,0,0,0,0,0) call total_energy(33,2,2,6,2,6,2,10,3,0,0,0,0,0,0,0,0,0,0) write(*,*)"As(33)" write(*,*)times write(*,*)total write(11,*)"As(33)" write(11,*)times write(11,*)total end subroutine Se34() use IMSL use global implicit none !---Se--- call sub_self(34,2,2,6,2,6,2,10,4,0,0,0,0,0,0,0,0,0,0) call total_energy(34,2,2,6,2,6,2,10,4,0,0,0,0,0,0,0,0,0,0) write(*,*)"Se(34)" write(*,*)times write(*,*)total write(11,*)"Se(34)" write(11,*)times write(11,*)total end subroutine Br35() use IMSL use global implicit none !---Br--- call sub_self(35,2,2,6,2,6,2,10,5,0,0,0,0,0,0,0,0,0,0) call total_energy(35,2,2,6,2,6,2,10,5,0,0,0,0,0,0,0,0,0,0) write(*,*)"Br(35)" write(*,*)times write(*,*)total write(11,*)"Br(35)" write(11,*)times write(11,*)total

end subroutine Kr36() use IMSL use global implicit none !---Kr--- call sub_self(36,2,2,6,2,6,2,10,6,0,0,0,0,0,0,0,0,0,0) call total_energy(36,2,2,6,2,6,2,10,6,0,0,0,0,0,0,0,0,0,0) write(*,*)"Kr(36)" write(*,*)times write(*,*)total write(11,*)"Kr(36)" write(11,*)times write(11,*)total end subroutine Rb37() use IMSL use global implicit none !---Rb--- call sub_self(37,2,2,6,2,6,2,10,6,1,0,0,0,0,0,0,0,0,0) call total_energy(37,2,2,6,2,6,2,10,6,1,0,0,0,0,0,0,0,0,0) write(*,*)"Rb(37)" write(*,*)times write(*,*)total write(11,*)"Rb(37)" write(11,*)times write(11,*)total end subroutine Sr38() use IMSL use global implicit none !---Sr--- 60

call sub_self(38,2,2,6,2,6,2,10,6,2,0,0,0,0,0,0,0,0,0) call total_energy(38,2,2,6,2,6,2,10,6,2,0,0,0,0,0,0,0,0,0) write(*,*)"Sr(38)" write(*,*)times write(*,*)total write(11,*)"Sr(38)" write(11,*)times write(11,*)total end subroutine Y39() use IMSL use global implicit none !---Y--- call sub_self(39,2,2,6,2,6,2,10,6,2,1,0,0,0,0,0,0,0,0) call total_energy(39,2,2,6,2,6,2,10,6,2,1,0,0,0,0,0,0,0,0) write(*,*)"Y(39)" write(*,*)times write(*,*)total write(11,*)"Y(39)" write(11,*)times write(11,*)total end subroutine Zr40() use IMSL use global implicit none !---Zr--- call sub_self(40,2,2,6,2,6,2,10,6,2,2,0,0,0,0,0,0,0,0) call total_energy(40,2,2,6,2,6,2,10,6,2,2,0,0,0,0,0,0,0,0) write(*,*)"Zr(40)" write(*,*)times write(*,*)total write(11,*)"Zr(40)" write(11,*)times

write(11,*)total end subroutine Nb41() use IMSL use global implicit none !---Nb--- call sub_self(41,2,2,6,2,6,2,10,6,1,4,0,0,0,0,0,0,0,0) call total_energy(41,2,2,6,2,6,2,10,6,1,4,0,0,0,0,0,0,0,0) write(*,*)"Nb(41)" write(*,*)times write(*,*)total write(11,*)"Nb(41)" write(11,*)times write(11,*)total end subroutine Mo42() use IMSL use global implicit none !---Mo--- call sub_self(42,2,2,6,2,6,2,10,6,1,5,0,0,0,0,0,0,0,0) call total_energy(42,2,2,6,2,6,2,10,6,1,5,0,0,0,0,0,0,0,0) write(*,*)"Mo(42)" write(*,*)times write(*,*)total write(11,*)"Mo(42)" write(11,*)times write(11,*)total end subroutine Tc43() use IMSL use global implicit none !---Tc--- 62