昨日行ったHe-He間距離に対する結合エネルギーカーブをRによってフィッティングしてみました。
6-12 LJ型の関数およびモース型関数に対してフィッティングを行いました。
- 6-12 LJ型関数
εは平衡距離での結合エネルギー(前回MO計算では 0.00152 kcal/mol)、r0は平衡距離(前回MO計算では 3.3856Å)
- モース型関数
Deは平衡距離での結合エネルギー、aはポテンシャル井戸の幅、reは平衡距離
CSV形式の入力ファイル
He間距離, エネルギーを以下のようなCSVファイル(C:/work/HeEne.csv)として用意する。
distance,energy 3.1,0.001686872 3.2,-0.000489081 3.3,-0.001352849 3.3856,-0.001519139 3.44,-0.001474021 3.5,-0.001350213 3.7,-0.000784827 4,-0.000227598 4.5,-1.43E-05 5,-4.39E-07 6,0
Rによるフィッティング実行
以下のRコマンドによりフィッティングを実行。
# CSVファイルを読み込む dat <- read.csv("C:/work/HeEne.csv") # データをプロット plot(dat, lwd=4) # 6-12 LJ関数をフィッティング用に定義 LJfrms <- as.formula(y ~ e*((r0/x)^12 - 2*(r0/x)^6)) # x, y のデータフレームに変換 ener <- data.frame(x=dat$distance, y=dat$energy) # フィッティング実行、start=list()にパラメータ初期値を与える # 初期パラメータが悪いと収束しない。その場合は初期パラメータを変えて再実行 Fit <- nls(LJfrms, start=list(e=0.00152, r0=3.3856), data=ener, trace=TRUE) # フィットされたパラメータを取得 # ()で囲むことによりパラメータがコンソールに表示される (e <- Fit$m$getPars()["e"]) (r0 <- Fit$m$getPars()["r0"]) # 3.0から6.0まで0.05刻みで出力する x <- seq(3.0, 6.0, 0.05) # 6-12 LJ関数を描画用に定義 LJfunc <- function(x, e, r0){e*((r0/x)^12 - 2*(r0/x)^6)} # フィットされたパラメータにより6-12 LJ関数を赤色でプロット、lwdは線の太さ lines(x, LJfunc(x, e, r0), col="red", lwd=2) # MO計算で求めたエネルギー・平衡距離によるプロット e <- 0.00152 r0 <- 3.3856 x <- seq(3.0, 6.0, 0.05) LJfunc <- function(x, e, r0){e*((r0/x)^12 - 2*(r0/x)^6)} # MO計算結果のパラメータにより6-12 LJ関数を青色でプロット lines(x, LJfunc(x, e, r0), col="blue", lwd=2) # モース関数をフィッティング用に定義 Morsefrms <- as.formula(y ~ D * (exp(-2*b*(x-re)) - 2*exp(-b*(x-re))) ) # フィッティング実行 Fit <- nls(Morsefrms, start=list(D=0.0015, b=3.0, re=3.3856), data=ener, trace=TRUE) # フィットされたパラメータを取得 (D <- Fit$m$getPars()["D"]) (b <- Fit$m$getPars()["b"]) (re <- Fit$m$getPars()["re"]) x <- seq(3.0, 6.0, 0.05) # モース関数を描画用に定義 MorseFunc <- function(x, D, b, re){D * (exp(-2*b*(x-re)) - 2*exp(-b*(x-re)))} # フィットされたパラメータによりモース関数を緑色の線でプロット lines(x, MorseFunc(x, D, b, re), col="green", lwd=4)