Self-consistency

0

P_pr(r)dr =

∫ _∞

0

r²R^∗_nlR_nldr = 1. (2.29) Therefore, we have the conclusion that if we want to normalize a radial wave function R_nl, we just compute the integral

∫ _∞

0

r²R^∗_nlR_nldr = k, k is a constant, and divided R_nl by√

k. Then the new function R⁰_nl= ^R√^nl

k is the radial wave function which has been normalized.

2.7 Self-consistency

Define M_{ef f} be the matrix which signifies the effective potential v_{ef f}, and M_{ef f} = M_ee+ M_xc+ M_r. (2.30) By (2.23), (2.24) and (2.30), the equation (1.18) can be represented as the form

(M_p+ M_l+ M_{ef f}) ~P_nl = e_nlP~_nl, where the vector ~P_nl is

13

P~_nl =











P_nl( ~r₁) P_nl( ~r₂) ...

P_nl( ~r₃₉₉)











.

It is an eigenvalue problem.

We can determine the correct electronic density and ~P_nl 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 H_kh (here H_kh = M_p + M_l + M_{ef f}). Then, computing the minimal eigenvalue emin of Hkh, and find an eigenvector ~P corresponding e_min. Normalizing ~P and we get the eigenvector ~P_out. Then we can compute the electronic density ~ρ_out. Now, if ~ρ_out = ~ρ_in it means that we get the right P_nl, otherwise, generate a new ~ρ_inas this ~ρ_out and construct a new H_kh. 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.

14

out 

^in

Table (2.7) YES

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

2. Compute ^_in, _{in i}| ^{in nl( )} |²

i occ

f P

 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( )} |²

i occ

f P

 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.

Figure (3.2)

Atomic number Idea data

–E(a.u.)

Data in this computation Errors(%)

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

18

Atomic number Idea data

-E(a.u.)

Data in this computation Errors(%)

36

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 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(%)

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

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.

20

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

40

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

44

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 sub_self(29,2,2,6,2,6,1,10,0,0,0,0,0,0,0,0,0,0,0)

在文檔中 原子基態能量的計算 (頁 20-99)