cms.integration.N_integration_metal のソースコード

import sys
import time
from math import exp, sqrt
import numpy as np
from scipy import integrate             # 数値積分関数 integrateを読み込む
from matplotlib import pyplot as plt


"""
金属の有限温度 T における電子密度 N の計算。
数値積分の制度の確認
"""


#定数
pi   = 3.14159265358979323846
h    = 6.6260755e-34    # Js";
hbar = 1.05459e-34      # "Js";
c    = 2.99792458e8     # m/s";
e    = 1.60218e-19      # C";
kB   = 1.380658e-23     # JK<sup>-1</sup>";


# Temperature
T = 300.0 # K

# Fermi energy
EF = 5.0 # eV

# EF近傍の nrange * kBTの範囲で積分
nrange = 6.0

# quad/romberg積分の相対積分精度、quad関数のGauss積分の最大次数/romberg積分の最大分割数
eps = 1.0e-8
maxorder = 100

#計算時間を計測する繰り返し数
ncycle = 300

# Treat argments
if __name__ == "__main__":
    argv = sys.argv
    if len(argv) >= 2:
        T   = float(argv[1])
    if len(argv) >= 3:
        EF = float(argv[2])

ekBT = e / (kB * T)    # eV単位のエネルギーを E/(kBT) に変換する係数
dE = nrange * kB * T / e

# 状態密度関数


def De(E):
    return sqrt(E)


# FE分布関数


def fe(E, T, EF):
    global kBT
    if T == 0.0:
        if E < EF:
            return 1.0
        else:
            return 0.0
    return 1.0 / (exp((E - EF) * ekBT) + 1.0)


# 被積分関数を定義


def Defe(E, T, EF):
    return De(E) * fe(E, T, EF)






def main():
    print("T = ", T, "k   EF = ", EF, " eV")
    print("")

    Ianaly = 2.0 / 3.0 * pow(EF, 1.5)
    print("Analytical integration for De(E) from ", 0, " to ", EF, " = ", Ianaly)
    ret = integrate.quad(De, 0.0, EF, epsrel = eps)
    print("quad   : ret=", ret)
    ret = integrate.romberg(De, 0.0, EF, rtol = eps, divmax = maxorder)
    print("romberg: ret=", ret)
    print("")

# 積分範囲
    (Emin, Emax) = (0.0, EF + dE)
    print("(1) Integration range: ", Emin, " - ", Emax, " eV")
    t0 = time.time()
    for i in range(ncycle):
        ret = integrate.quad(lambda E: Defe(E, T, EF), Emin, Emax, epsrel = eps)
    dt = time.time() - t0
    print("quad   : ret=", ret, "  time=", dt, " s")
    t0 = time.time()
    for i in range(ncycle):
        ret = integrate.romberg(lambda E: Defe(E, T, EF), Emin, Emax, rtol = eps, divmax = maxorder)
    dt = time.time() - t0
    print("romberg: ret=", ret, "  time=", dt, " s")
    print("")

# 積分範囲
    (Emin, Emax) = (0.0, EF - dE)
    print("(2) Integration range: ", Emin, " - ", Emax, " eV")
    Ianaly = 2.0 / 3.0 * (pow(Emax, 1.5) - pow(Emin, 1.5))
    print("Analytical integration for De(E) from ", Emin, " to ", Emax, " = ", Ianaly)
    t0 = time.time()
    for i in range(ncycle):
        ret = integrate.quad(lambda E: Defe(E, T, EF), Emin, Emax, epsrel = eps)
    dt = time.time() - t0
    print("quad   : ret=", ret, "  time=", dt, " s")
    t0 = time.time()
    for i in range(ncycle):
        ret = integrate.romberg(lambda E: Defe(E, T, EF), Emin, Emax, rtol = eps, divmax = maxorder)
    dt = time.time() - t0
    print("romberg: ret=", ret, "  time=", dt, " s")
    print("")


# 積分範囲
    (Emin, Emax) = (EF - dE, EF + dE)
    print("(3) Integration range: ", Emin, " - ", Emax, " eV")
    t0 = time.time()
    for i in range(ncycle):
        ret = integrate.quad(lambda E: Defe(E, T, EF), Emin, Emax, epsrel = eps)
    dt = time.time() - t0
    print("quad   : ret=", ret, "  time=", dt, " s")
    t0 = time.time()
    for i in range(ncycle):
        ret = integrate.romberg(lambda E: Defe(E, T, EF), Emin, Emax, rtol = eps, divmax = maxorder)
    dt = time.time() - t0
    print("romberg: ret=", ret, "  time=", dt, " s")
    print("")



if __name__ == '__main__':
    main()