• 沒有找到結果。

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

N/A
N/A
Protected

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

Copied!
174
0
0

(1)

## 學

2

2

(2)

2

2

(3)

(4)

2

2

(5)

(6)

kin

ee

ec

xc

(7)

2

2

(8)

2

(9)

2

2

(10)

G

G

ec

m

H

ec

ext

ext

m

H

m

H

0

ee

xc

G 0

(11)

i

pl

lm

2

i i

i occupied

### density of O-atom is

2 2 2 2 2 2

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

G

G 0

ee

ec

xc

ee

ee

ee

(12)

2 * * 0

### 2

i i i i i i i occupied i occupied

∈ ∈

ec

ec

ext

xc

xc

xc

xc

xc

,

I J cc I J I J

(13)

1 1 2 2 n n i i

i

i

i

i

(14)

2

2

1 2

n

(15)

pl

2

(1,0,0)

(2,0,0)

(2,1,-1)

(2,1,0)

(2,1,1)

2

2

2

2

(16)

1 1

i

i i j i j j

− =

1

2

10

i i

i

2 2 * ,

i i

i j

i j i i j

2

2 * ,

i j i j i j

* * ,

i j i j i j

* * * *

i

i

i

i

(17)

j

j

j

j

* * * *

i j i j i j i

j

* *

*

i i

i

j j

j

*

*

*

i i

i

j j

j

* *

*

i i

i

j j

j

3

ee

ee

* ,

i j ee i j i j

(18)

2 2

i i k k i k ee

* * , , ,

### '

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

*

*

i j k h i j k h

i j

-I ext I I

I

I

2

ec ext i i i

* , i j ext

i j i j

(19)

I

0 0

0 0

N N

i

ext

j

* * [ , ]\A-

### ( ) ( ) ( )

ext i j ext i j a b part part

(20)

xc

xc

* ,

i j xc i j i j

xc

xc

x

c

xc

x

c

1 3

s

x

s

1 2

s s c s s s s

1

2

(21)

* ,

i

j i j i j

1 2

10

* ,

i j i j i j

12

22

102

1 2 10

12

22

102

1 2

10

12

22

102

### Therefore, we have

(22)

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

1 2

10

kin

ec

ee

xc

kin

ec

ee

xc

i i

10 10 2 2 2 1 1

### 2

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

=

=

10 10 1 1

### 2

ec

j i ext j j ext i j i ext j

j j i

=

=

i j k h k h i j

### 

(23)

10 10 , , 1 , , 1 10 10 , , 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

= = = =

### 

10 10 , , 1 , , 1

### '

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

= =

10 , , 1

### '

i j k i j k h i j k

=

1 2

10

i

ee

j

10 1

j i ee j

j

=

10 10 1 1

### 2

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

=

=

10 2 1

j i

j

### +2

i ext j i ee j i xc j

j i

=

10 i j 1

j

j i

=

### ∑

(24)

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

j j j j j j

= = =

### 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

old

2 2 2

1 2 10

new

1 2 10

new

old

New Old New

(25)

cc

old

new

(26)
(27)
(28)
(29)

2222

(30)
(31)

(32)
(33)

(34)
(35)
(36)
(37)

2

2

2

(38)

*

2

ext

ee

xc

i

i

i

i

i

c

c

c

c

(39)

### 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

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

do h = 1,basics

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

GCG method is developed to minimize the residual of the linear equation under some special functional.. Therefore, the minimality property does not hold... , k) in order to construct

Understanding and inferring information, ideas, feelings and opinions in a range of texts with some degree of complexity, using and integrating a small range of reading

Writing texts to convey information, ideas, personal experiences and opinions on familiar topics with elaboration. Writing texts to convey information, ideas, personal

Writing texts to convey simple information, ideas, personal experiences and opinions on familiar topics with some elaboration. Writing texts to convey information, ideas,

Wang, Solving pseudomonotone variational inequalities and pseudocon- vex optimization problems using the projection neural network, IEEE Transactions on Neural Networks 17

3: Calculated ratio of dynamic structure factor S(k, ω) to static structure factor S(k) for &#34;-Ge at T = 1250K for several values of k, plotted as a function of ω, calculated

Define instead the imaginary.. potential, magnetic field, lattice…) Dirac-BdG Hamiltonian:. with small, and matrix

Total energies and Cartesian coordinates for the lowest singlet and triplet states of n-acenes (n = 2 to 46) by spin-unrestricted TAO-LDA (θ = 7 mHartree)/6-31G ∗ (S3 to S88)..