有機超伝導体や非線形光学材料として分子性結晶に寄せられる期待は大きい。これらの材料としての可能性は結晶の構造から評価できるが、理論的な予測法が確立されれば、新規の分子性結晶の開発に有用な情報を与えることができる。

　分子性結晶の構造の予測には次の2点がポイントとなる。

　1)分子間の相互作用を適切に記述するポテンシャル関数の決定。

　2)与えらたポテンシャル関数の下での最も安定なパッキングの探索。

　未知分子の結晶構造を予測することを考えると、この2つのプロセスが実験的データを参照せずに行えることが望ましい。この観点から、本研究では、まず2)に関して、空間群や格子定数等の実験的情報を必要としない、パッキングの予測法の開発を行った。次に1)に関して、分子軌道計算に基づいて決定したポテンシャル関数を用いて結晶構造予測を行い、分子式のみから結晶構造を予測できることを示した。また、結晶構造の解析に実験的な困難を伴う場合には、理論的な構造予測法が支援となる。この例として、高圧下のCO₂分子の結晶構造と、低温で結晶化する1,2ジメトキシエタン分子の結晶構造の予測を行った。

　<結晶構造の探索> 分子性結晶の構造の予測は、分子が最も安定に集合した構造を探索することを方針とする。構造の安定性の熱力学的な指標は自由エネルギーであるが、ここではポテンシャルエネルギーで近似し、これを最小にするような構造(グローバルミニマム)を周期的境界条件を課した8個の分子より成る系の中から探索する。また、未知分子の結晶構造を予測する立場から、しばしば用いられる空間対称性の制約は、探索対象の一般性を保つために設けない。

　<アルゴリズム> グローバルミニマムを得るためには網羅的な探索法が必要であるが、このような小さな系でもポテンシャルの極小点は非常に多いことが計算上の難点である。そこで、グローバルミニマムを効率よく求める為に、分子動力学(MD)を用いた以下のような方法を考案した。

　まず、ポテンシャル曲面のエネルギーの低い領域を走査するため、ランダムな初期構造から出発し、特定のローカルミニマムには捕われない程度の低温で定圧MD計算を行う(シェーキング)。その軌跡上でのポテンシャルエネルギーの極小点は、パッキングの良い構造としてサンプリングされる。この際、ポテンシャル局面の同一の谷に属する構造を重複して採らないように、原子間距離の標準偏差によって類似性を調べながらサンプリングに間隔を置く。次に、サンプリングされた構造を初期値としたエネルギー極小化計算によって、ポテンシャル極小点を求める(クエンチング)。得られた極小点の構造を、エネルギー、規約格子へ変換した時の格子定数、空間対称性に関して分類する。ポテンシャル関数の不正確さを考慮して、独立な極小点の構造のうちエネルギーの低い方から数個を結晶構造の候補とする。この方法は、プログラムMDCP(Molecular Dynamcs for Crsytal Packing)にコーディングした。

　<テスト計算> CO₂、ベンゼン、ピリミジン分子を対象としてテスト計算を行った。分子は結晶中の構造を持つ剛体として近似し、分子間のポテンシャル関数には、Williams等のものを用いた。

　CO₂,ベンゼン分子では、200Kでの2万ステップのシェーキングからそれぞれ98,40個のサンプル構造が得られ、更にクエンチングと分類の後には約10個の異なる候補構造となった(図1)。この中で最もポテンシャルエネルギーの低い構造は、実測の結晶構造とほぼ等しいものであった(表1)。本研究の方法で、グローバルミニマムの探索が効率よく行われ、結果として空間対称性や格子定数も含めての結晶構造の予測が可能となることが示された。

図表表1.結晶の候補構造の格子定数 / 図1.クエンチ後の構造のポテンシャルギーの分布(右)

図表

　一方ピリミジン分子では、10万ステップのシェーキングからも結晶構造に対応する極小点は得られていない。これは、グローバルニマムとポテンシャルエネルギー値の近い極小点が多く存在し、これらを探索しきれなかったことが原因である。このような分子には多大な計算時間が必要となる。

　実験的データを参照せずにポテンシャル関数を決定する際には、分子配置毎に相互作用エネルギーを求めることのできる量子化学的手法が有力である。ここではCO₂分子を例にとり、ab initio分子軌道計算に基づいてポテンシャル関数を決定し、完全に非経験的な結晶構造予測を試みる。

　<ポテンシャル関数の決定> CO₂分子の変形は考えないこととし、MP2/6-311G^*レベルの計算より得られる分子構造に固定した。分子間相互作用エネルギーは、atom atom型のポテンシャル関数で表わし、各原子間のポテンシャル関数型には、exp-6-1型を用いた。-1乗項のパラメータである部分電荷は、孤立分子に対するMP2/6-311G(2d)レベルの計算での静電ポテンシャルにフィットした。CO₂2量体の異なる40配置について相互作用エネルギーをMP2/6-311G(3d)レベルで計算し、これを再現するようexp-6項のパラメータを最小2乗法で決定した。

　<ポテンシャル関数の精度> 上のポテンシャル関数を用いて、結晶の格子エネルギーと構造の計算を行い、精度を評価してみた(表2)。実測値から求めたKuchtaらのポテンシャル関数には劣るが、格子エネルギーや、低圧相(0GPa)と高圧相(20GPa)の格子定数は、数%の誤差で実測値を再現する。ab initio計算で求めた2分子間の相互作用の総和で全相互作用が適切に記述され、全エネルギーの規模や極小点の位置が正しく与えられることが分かる。

表2.格子エネルギーと格子定数の計算値

　<結晶構造予測>上のポテンシャル関数を用い、低圧(0GPa)及び高圧(20GPa)のCO₂の結晶構造について、1.の手法によってパッキングの予測計算を行った。両圧力下で最低のエンタルピーを示したのは、実測構造に対応するエネルギーミニマムであった。ab initio分子軌道計算から求めたポテンシャル関数を1.で述べた手法に用いることによって、完全に非経験的な結晶構造の予測が可能であることが示された。

　CO₂の結晶は、常圧ではドライアイスとして知らる立方晶(Pa3)であるが、12GPa付近まで加圧すると他の構造へ相転移を起こす。高圧相の構造について、最近の粉末X線回折の実験では、斜方晶(Cmca)を支持する回折パターンが得られている。この構造への転移について、熱力学的安定性と動力学的安定性の観点から考察を行った。

　<自由エネルギー計算> CO₂分子は剛体として近似し、分子間の相互作用には2.で求めたポテンシャル関数を用いた。1.の方法を応用し、0から20GPaの外圧下でのエンタルピーミニマムを求め、その自由エネルギーを調和近似を用いて計算した。低圧では立方晶(Pa3)のエネルギーが最も低いが9GPaからは斜方晶(Cmca)のエネルギーが最低となった。斜方晶(Cmca)への転移が熱力学的な安定性の面から説明できた。

　<MDシミュレーション> 108分子の系に対する分子動力学計算を行ったところ、実験よりも高温ではあるが1000K,20GPaで斜方晶(Cmca)への転移がみられた。斜方晶(Cmca)への転移は低圧相(Pa3)から最も容易なものとしても理解できる。また、等方的な12GPaの外圧の他に1GPaの(1,1)成分の外部応力を付加すると、実験の温度(295K)でも斜方晶(Cmca)への転移がみられる。圧力印加の際の、試料内の圧力の不均一が転移の引き金となっていると考えられる。

　1,2ジメトキシエタン(DME)は、気相や液相中では種々の配座を取って混在しているのに対し、結晶中ではTGT型しか取らない。結晶構造に関するX線解析の報告はないため、ここでは、1.で開発した方法を応用してその結晶構造を調べてみた。

　<計算方法> 計算に用いる分子は剛体として近似し、その構造には孤立状態で配座エネルギーの低いTTT,TGT,TGG’型を選んだ。ローカルミニマムが多いことが予想されるため、分子性結晶に出現頻度の大きい9種の空間対称性を初期構造に課して探索範囲を絞った。これより出発し、100Kのシェーキングを1000ステップ行うことで空間対称性をある程度崩した後、クエンチングを行って候補構造を作成した。

　<ポテンシャル関数> MM3のポテンシャル関数をVDW相互作用に用い、更に非結合相互作用を改善するため、分子軌道計算より得られる静電ポテンシャルにフィッティングした部分電荷を各原子上に置いて静電相互作用を取り込んだ。

　<計算結果> 3つの配座異性体から得られる構造の分子間ポテンシャルエネルギーの最小値を比較すると、TGT型のものが最も低い。分子内のポテンシャルエネルギーを考慮して全エネルギーについて比較しても、TGT型のものがTTT,TGG’型のものに対してそれぞれ0.92,3.94kJ/mol低い。したがって、分子が集合することよって固体中ではTGT型が有利になることが、計算からも裏付けられたことになる。この構造の空間対称性はPnaaで、分子は互いの双極子の方向を揃えてクーロン相互作用に都合の良いパッキングをしていて、もっともらしい結晶構造となっている。