量子力学はLennard-Jonesポテンシャルを再現できるか
原理的にはできないと困るわけだけど、化学的精度と言われる1kcal/mol程度の計算は現在でも多くの場合容易ではない。でも、系が簡単ならできるのではという期待も微妙にある。単純に、二分子系で構造最適化して、ポテンシャルの深さと平衡分子間距離からLJパラメータを決めても、無極性分子同士の分子間相互作用は弱いので、構造最適化の精度がでない可能性は高い。その場合は、分子間距離を変えて一点計算を繰り返すほうがよいかもしれない。まぁ、既に何人もの人がやってるだろうけど。計算には、Gamessを使用
LJポテンシャルは
で、ポテンシャルが極小になるrを$r_m$とすると、極小値は-εで
という関係がある。
まず、一番簡単なヘリウム。インプットは以下のように適当に作った(とりあえず、HFで、6-31G(d)基底)
! rungms He2.inp 11-64 1 0 > He2.log $contrl runtyp=optimize scftyp=rhf coord=unique $end $basis gbasis=n31 ngauss=6 ndfunc=1 $end $guess guess=Huckel $end $statpt OPTTOL=1.0e-9 HSSEND=.T. $end $data helium dimer C1 He 2 0.0 0.0 0.0 He 2 0.0 0.0 2.2 $end
構造最適化の結果、分子間距離は3.2345263Åで、エネルギーは-5.7103221681(hartree)となる。σ=2.88Å程度。ヘリウム一原子のエネルギーは同条件の計算で、-2.8551604262(hartree)らしいので、ε=1.3157000005037389e-6(hartree)=0.00083(kcal/mol)。ヘリウムのLJパラメータは、見る文献によって、意外とバラツキがあって、どれがよいのか分からないけど、σが2.6Å前後で、εが0.02kcal/mol程度なのは確かっぽいので、距離はともかく、エネルギーが合ってない
同じく、HeでCI計算をやってみる
$contrl runtyp=optimize scftyp=rhf CITYP=GUGA $end $CIDRT NFZC=0 NDOC=2 NVAL=8 IEXCIT=2 $END $GUGDIA PRTTOL=0.00001 $END $basis gbasis=n31 ngauss=6 ndfunc=1 npfunc=1 $end $statpt OPTTOL=1.0e-9 $end $data helium dimer C1 He 2 0.0 0.0 0.0 He 2 0.0 0.0 3.3 $end
r_m=3.0420800Åが平衡距離で、σ=2.71Åとなって、距離は大分近づいた。この時-5.7742572988(hartree)で、一原子のエネルギー-2.8873650277(hartree)なので、だめっぽい。仕方ないので、100Å離して一点計算すると-5.7742486047(hartree)らしい。なので、ε=0.0055(kcal/mol)
HeでCCSD(T)
$contrl runtyp=optimize scftyp=rhf CCTYP=CCSD(T) NUMGRD=.T. $end $basis gbasis=n31 ngauss=6 ndfunc=1 npfunc=1 $end $statpt OPTTOL=1.0e-9 $end $data helium dimer C1 He 2 0.0 0.0 0.0 He 2 0.0 0.0 3.3 $end
r_m=3.0416303Åが平衡距離で、σ=2.71Å。エネルギーは-5.7747387998(hartree)で、一原子のエネルギー-2.8873650277(hartree)なので、ε=8.7443999996494881e-06(hartree)=0.0055(kcal/mol)。CI計算とほぼ同一の結果
cc-pVDZ基底を持ってきて、CCSD(T)計算
! cc-pVDZ basis ! from https://bse.pnl.gov/bse/portal $contrl runtyp=optimize scftyp=rhf CCTYP=CCSD(T) NUMGRD=.T. $end $statpt OPTTOL=1.0e-9 $end $data helium dimer C1 He 2 0.0 0.0 0.0 S 3 1 38.3600000 0.0238090 2 5.7700000 0.1548910 3 1.2400000 0.4699870 S 1 1 0.2976000 1.0000000 P 1 1 1.2750000 1.0000000 He 2 0.0 0.0 3.3 S 3 1 38.3600000 0.0238090 2 5.7700000 0.1548910 3 1.2400000 0.4699870 S 1 1 0.2976000 1.0000000 P 1 1 1.2750000 1.0000000 $end
r_m=3.0907625Åで、σ=2.75Å。エネルギーは-5.7751959267(hartree)で、50Å離した時の一点計算の結果-5.7751896622(hartree)を基準とすると、ε=0.004(kcal/mol)となる。まぁ、流石に厳しいのか
Arでもやってみる。ただのHF
! rungms Ar2.inp 11-64 1 0 > Ar2.log $contrl runtyp=optimize scftyp=rhf coord=unique $end $basis gbasis=n31 ngauss=6 ndfunc=1 $end $guess guess=Huckel $end $statpt OPTTOL=1.0e-9 HSSEND=.T. $end $data argon dimer C1 Ar 18 0.0 0.0 0.0 Ar 18 0.0 0.0 3.2 $end
分子間距離は4.5129590Åで、σ=4.02Å。エネルギーは-1053.5475039624(hartree)。一原子のエネルギーは-526.7737449209(hartree)で、ε=1.412e-5(hartree)=0.00836(kcal/mol)。実測値は、σ=3.4Å,ε=0.2kcal/molなので、微妙
MP2計算にすると(基底関数は同じ)、分子間距離は4.1782062Åで、σ=3.722Å。エネルギーは-1053.8222010388(hartree)。一原子のエネルギーは、-526.9110529280(hartree)より、ε=0.06(kcal/mol)で、改善される
折角なので、二原子分子N2でもやってみる。まず、一分子の構造最適化をする
$contrl runtyp=optimize scftyp=rhf mplevl=2 $end $basis gbasis=n311 ngauss=6 ndfunc=1 $end $guess guess=Huckel $end $statpt OPTTOL=1.0e-9 HSSEND=.T. $end $data nitrogen molecule C1 N 7.0 0.0000000000 0.0000000000 0 N 7.0 0.0000000000 0.0000000000 1.2 $end
トータルエネルギーは-109.3027176037(hartree)。二原子分子の場合、相対的な配置も問題になるけど、とりあえず平行に置いてみる。でまぁ、適当にインプットを作る。この程度の計算なら、いきなり重めの基底関数から始めても余裕
$contrl runtyp=optimize scftyp=rhf mplevl=2 $end $basis gbasis=n311 ngauss=6 ndfunc=1 $end $guess guess=Huckel $end $statpt OPTTOL=1.0e-9 HSSEND=.T. $end $data N2 dimer C1 N 7.0 0.0000000000 0.0000000000 0.0398482444 N 7.0 0.0000000000 0.0000000000 1.1601517556 N 7.0 3.0000000000 0.0000000000 0.0398482444 N 7.0 3.0000000000 0.0000000000 1.1601517556 $end
これで最適化して、基準振動を見ると、虚のモードが出てしまう。というわけで、別の配置からやり直し
$contrl runtyp=optimize scftyp=rhf mplevl=2 $end $basis gbasis=n311 ngauss=6 ndfunc=1 $end $guess guess=Huckel $end $statpt OPTTOL=1.0e-9 HSSEND=.T. $end $data N2 dimer C1 N 7.0 0.0000000000 0.0397494066 0.0000000000 N 7.0 0.0000000000 1.1602505275 0.0000000000 N 7.0 3.0000000000 0.0000000000 0.0397494724 N 7.0 3.0000000000 0.0000000000 1.1602505935 $end
これだと、虚のモードは出ない。最適化後の配置は
N 7.0 -0.1514279733 0.1907019398 0.0836249847 N 7.0 -0.5073563096 1.1933370214 -0.2673972831 N 7.0 3.1529363171 0.0835077242 0.1898945065 N 7.0 3.5058479658 -0.2675467513 1.1938778579
分子の重心間距離で見ればr_m=3.82Åで、σ=3.40Å。エネルギーは-218.6060802776(hartree)なので、ε=0.40(kcal/mol)程度。
http://www.nucleng.kyoto-u.ac.jp/people/ikuji/edu/ebeam/lennard.html
によると、σ=3.68Å,ε=0.18kcal/molで、まぁぼちぼち
計算方法が行き当たりばったりで統一性がないけど、まぁ、ちゃんとやればそれなりの結果が出そう？
課題:
・相互作用エネルギーが小さすぎて、BSSEなどの数値誤差の影響がでかい
・実験値は、そもそも二体ポテンシャル近似とか過程を置いた上でLJパラメータを出しているので、純粋な二分子系の計算結果とどれくらい一致するのか不明
・多原子分子の場合、LJポテンシャルは、距離のみの関数なので、ポテンシャルエネルギーは、最適化した配置ではなく、(分子の重心を固定した上で)分子同士の相対配置全部に渡って平均化したエネルギーを考えるのが正しい気もするし、そもそもエネルギー極小となる配置が"本質的に複数"(つまり、一様な回転や並進操作によって互いに移らない)存在する可能性もある
