読者です 読者をやめる 読者になる 読者になる

生物物理計算化学者の雛

主に科学に関する諸々を書き留めています。

Rにより結合エネルギー関数フィッティングを行う

R 分子シミュレーション

昨日行った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)
フィッティング結果


6-12 LJ関数でフィッティングを行うことで、近距離(<3.4Å)でのポテンシャルの立ち上がりが遅い部分は改善されました。しかしながらMO計算結果ではHe間距離を長くするとすぐにエネルギーが0に収束する一方、6-12 LJ関数では距離を長くしてもなかなか収束しない部分はさほど変化がありません。
それに対してモース関数では近距離・遠距離両方で見事にMO計算結果を再現できていることがわかります。