JP5011689B2 - 分子シミュレーション方法及び装置 - Google Patents

Description

本発明は、分子動力学シミュレーションあるいは分子モンテカルロ法の手法により分子シミュレーションを行うための方法及び装置に関し、特に、非結合相互作用の高速計算に用いるマルチグリッド法の改良に関する。

分子動力学シミュレーションの全般的な技術：
生体分子の分子動力学（ＭＤ；Molecular Dynamics）シミュレーションは、分子シミュレーションとも呼ばれ、疾病の治療薬や有用酵素などを分子設計するに使用されている重要な方法論である。この分子シミュレーションで用いられているモデルは、いわゆる分子力学（ＭＭ；Molecular Mechanics）モデルである。この分子力学モデルを簡単に説明すると、対象とする分子であるタンパク質や水分子を、バネで繋がれたビーズ（珠）で表現し、ビーズすなわち分子を構成する原子には電荷が割り振られているものである。分子力学モデルを用いる場合、分子系のエネルギーは、結合相互作用と非結合相互作用に大別される。このうち、結合相互作用は、結合長・結合角・結合捩れ角などの変化に伴う相互作用である。一方、非結合相互作用には、クーロン(Coulomb)相互作用やファン・デル・ワールス(van der Waals)相互作用が含まれる。

分子シミュレーションでは、分子における初期配置を定め、分子を構成する各原子に電荷を割り振って初期状態とし、上述した結合相互作用及び非結合相互作用を介して各分子が初期状態からどのように運動（例えば変形）し、またそれに伴って系のエネルギーがどのように変化していくかを計算する。多数の初期配置から分子シミュレーションを実行することによって、最終的に最も安定な分子の配置を定めることができ、これは、例えば、タンパク質とそのタンパク質に対してリガンドとして作用する薬剤との相互作用を推定したり、タンパク質自体の機能を推定する際の重要な情報となる。

典型的なタンパク質分子系の分子シミュレーションを行う場合に、大部分の計算時間を占めるのが非結合相互作用、特に、クーロン相互作用の計算である。考察対象系に含まれる原子に順番に番号をふり、原子の帯びている電荷（部分電荷）をｑ_iと表すことにすると、クーロン相互作用エネルギーＵは次式で表される：
Ｕ＝Σ_i Σ_j ｑ_iｑ_j／ｒ_ij (1)
ここで、ｒ_ijは原子ｉと原子ｊの間の距離であり、ｉ，ｊに関する和計算（サメイション）は、ｉ＜ｊの条件下で行うものとする。部分電荷の値は、分子軌道（ＭＯ；Molecular Orbital）法などの理論的手法や、実験データとのフィティングに基づく経験的手法で決定される。系の原子数をＮとすれば、クーロン相互作用の計算量は、オーダーでＮ²回となる。

クーロン相互作用が特に重要なのは、クーロン相互作用が距離の２乗に反比例する長距離型の相互作用であって、相互に作用しあう粒子ペアの個数が莫大なものとなるからである。これに対してファン・デル・ワールス相互作用は、一般に距離の６乗に反比例する短距離型の相互作用であるから、計算対象の粒子数がずっと少なくなり、その計算量もクーロン相互作用の場合に比べてずっと小さくなる。

前述の分子シミュレーションを行って得られた結果から有意義な結果を抽出するには、全分子にわたる計算を１回行うことを１ステップとして、少なくとも１００万ステップ以上、望むらくはその千倍以上の計算を繰り返す必要があり、非常に莫大な計算時間が必要である。この理由のために、薬の設計などで分子シミュレーションを活用するために、非結合相互作用の高精度・高速計算技術には大きな要請があった。なお、結合相互作用は、化学結合を介して直接つながっている原子間の相互作用であるため、非結合相互作用に比べて計算対象となる粒子ペア数は極めて少なくなる。

このように非結合相互作用に関して莫大な計算をしなければならないという要請に応えるために、分子シミュレーションのためのソフトウェアに関しては、様々なアルゴリズムが提案されて実装されている。また、分子シミュレーションを実行するためのハードウェアとしても、いわゆるＰＣクラスタや、地球シミュレータ、あるいはＩＢＭ社製ＢｌｕｅＧｅｎｅが活用されている。専用ハードウェアとしては、富士ゼロックス社製のＭＤエンジンや理化学研究所のＭＤＧｒａｐｅ等が開発されている。しかしながら、ソフトウェア手法、ハードウェア手法のいずれをとっても、分子シミュレーションの実行の観点からはいまだに高速化の度合いは不十分であり、さらなる高速化に対する要請が高い。

分子シミュレーションにおける代表的な２種類の境界条件：
分子シミュレーションで頻繁に使用される境界条件として代表的なものは、周期境界条件と真空境界条件の２種類がある。

周期境界条件は、例えば結晶のように、対象とする分子が周期的に配置しているというモデルに対応するものである。周期境界条件では、周期的に配置しているということからＦＦＴ（高速フーリエ変換）などの計算手法を使用でき、これによって計算速度を向上させることができる、とされている。

真空境界条件はいわゆる孤立系に対する計算であり、対象となる分子（例えばタンパク分子）が真空中に置かれた大きな球状の溶媒（典型的には水）の中に浸されている系がその典型的な例として挙げられる。真空境界条件は、一般には計算速度が遅いとされているものの、分子におけるコンフォメーション変化などを迅速にシミュレーションできるとも言われている。

周期境界条件、真空境界条件のいずれも人為的なものであって、それぞれ適不適があり、また課題を有することが指摘されており、そのため、周期系（周期境界条件）と孤立系（真空境界条件）の両方の境界条件に適用できる計算アルゴリズムが要望されている。

周期境界条件に対する既存アルゴリズムと問題点：
周期境界条件に対する効率的なアルゴリズムとしては、Dardenら(1995)のＰＭＥ(Particle Mesh Ewald)法という方法（非特許文献１参照）が著名である。この方法は、ポテンシャルを、粒子間で直接計算する部分(Particle-Particle)）と残りの滑らかな部分とに分離し、滑らかな部分に関しては、粒子をメッシュ上にマッピングし、畳み込み(convolution)型の計算を３Ｄ−ＦＦＴ（３次元高速フーリエ変換）を用いて高速計算することによって、最終的に粒子−メッシュ間(Particle-Mesh)相互作用を、効率よく計算する方法である。

このＰＭＥ法は、バイオシミュレーション分野で広く流布しているＭＤコード（分子シミュレーション用計算機コード）“Ａｍｂｅｒ”はもとより、他の様々なＭＤソフトでも実装されている優れた方法である。１６〜３２台程度までのＰＣクラスタでは、ＰＭＥ法のパラレル化性能（演算の並列処理化性能）は非常に良好である。しかし、それ以上の大規模になるとパラレル化率は低下し有効ではなくなるという問題も残っている。その原因は、大量のＣＰＵに対しては３Ｄ−ＦＦＴ演算のパラレル化が難しく、パラレル化率が低下してしまうためである。このパラレル化率低下を回避するために、ＩＢＭ社のＢｌｕｅＧｅｎｅでは３Ｄ−ＦＦＴを高速計算するためのハードウエア設計上の特別な配慮がなされている。

真空境界条件に対する既存アルゴリズムと問題点：
一方、真空境界条件（すなわち孤立系）に適用できる効率的なアルゴリズムとしては、Greengardらの高速多重極子展開法（Fast Multipole Method法、ＦＭＭ法）があり（非特許文献２参照）、“Ａｍｂｅｒ”とは別の著名なＭＤコードである“Ｃｈａｒｍｍ”などで実装されている。しかし、ＦＭＭ法は、非常に多数の高次多重極子の展開項を計算せざるをえないために、計算量がオーダーＮの計算手法とはいえ、計算がスローになっている（Ｎは系に含まれる原子の個数）。並列計算に用いるＣＰＵ台数が１６〜３２以下であっても、孤立系では３Ｄ−ＦＦＴが使えないために、むしろ周期系よりも取扱いが難しく、高速化しにくくなっている。実際、ＭＤコード“Ａｍｂｅｒ”では高速多重極子展開法は実装されておらず、多くのユーザーは、２つの粒子間の距離が一定距離以上離れたものについてはそれらの２粒子間の非結合作用を無視（カットオフ）して計算を行うか、もしくはＭＤエンジンなどの専用ボードを使用してカットオフなしの厳密計算を行っている。

なお、水などの球状の溶媒とその中のタンパク質からなる系のように孤立した分子系に対して周期系の演算方法を適用できるようにするために、そのような孤立分子系を周期的に配置した系を考えることもできる。周期的に配置していれば、ＦＦＴなどの適用も可能になる。しかしながらその場合、本来孤立しているはずの分子が周期的に存在していることによって生じる相互作用を回避するために、従来は、周期構造における分子と分子との間隔を、分子系の大きさよりも大きく設定する必要があり、結果として、より多点数のＦＦＴを実行しなければならず、計算量が大きくなるという問題が生じる。そのため、このように孤立分子系を周期的に配置することのアイディアは提案されているものの（非特許文献６）、実装されている例は稀である。

ＭＤ専用ボードの登場と問題点：
ＭＤ計算における粒子−粒子の相互作用を直接的に高速計算する専用ハードウェアとして、各種のＭＤ専用ボードが開発・市販されている。前述のＰＭＥやＦＭＭなどのアルゴリズム面で工夫された方法論と専用ボードとを組み合わせた計算もされているが、ＰＭＥやＦＭＭには専用ボードで計算不能な部分の割合が多く残っているので、このアルゴリズムを併用する場合には、併用しない場合に比べて高速化の程度が低下してしまう。その結果、ＰＭＥ法やＦＭＭ法以外の、専用ボードに適した手法に対する強い要請が出てきている。

マルチグリッド法の登場と問題点：
このような高速化・大規模パラレル化に対する強い要請に基づいて、最近、Skeelらは多重格子法（マルチグリッド(Multigrid)法）に基づく高速化を提案した（非特許文献３参照）。マルチグリッド法は、分子力学モデルにおけるファン・デル・ワールス・パラメータに対して特殊な混合則を採用する場合には、クーロン力に加えてファン・デル・ワールス力の計算にも適用可能となる。しかし後述するように、Skeelらの方法は高速ではあるが低精度という大きな弱点があり、実用上、十分な精度を得ることはできなかった。

以下、Skeelらのアルゴリズムを３項目に分けて説明する。

ポテンシャルを長距離型と近距離型の部分に分離：
Skeelらのアルゴリズムの基本的アイデアは１９７０年台より開発されてきたマルチグリッド法を応用したものである。本来計算したいポテンシャル関数が粒子間距離ｒの関数として、クーロン型（すなわち距離ｒに反比例）である場合で説明すると、図１ａに例示されているように、もとのポテンシャル関数Ｇ（ｒ）（＝１／ｒ）を、原点において特異性がない滑らかな関数Ｇ_R（ｒ）と、その残りを表す関数Ｇ^S _R（ｒ）との２つに分離する：
Ｇ（ｒ）＝Ｇ_R（ｒ） ＋ Ｇ^S _R（ｒ）
＝（長距離型）＋（短距離型） (2)
ここでＲは、遮蔽半径と呼ばれるパラメータである。長距離型関数Ｇ_R（ｒ）に必要な性質は、ｒが大きくなるにつれて元の関数Ｇ（ｒ）＝１／ｒに滑らかに収束し、特にｒが距離Ｒ以上では元の関数に完全に一致し、さらに、原点近傍ではＲが大きいほど平滑化されているという性質である。その結果、短距離型関数Ｇ^S _R＝Ｇ−Ｇ_Rは、ｒが距離Ｒ以下のときのみ非ゼロ値をとり、距離Ｒにおいて滑らかにゼロに収束することになる（図１ａ参照）。短距離型関数は粒子−粒子間相互作用として扱われ、長距離型関数は後述する格子を経由して扱われることになる。上添字のＳは、Short-Range（短距離）およびSplit（分離）を示している。このポテンシャル関数を、以下では慣例にしたがってカーネルと呼ぶ。

階層的な格子構造：
一般に、分子シミュレーションでは、考察粒子系を格子（あるいは、細胞、メッシュ、グリッド、セルなどと呼ばれる）で分割して、格子ごとに処理を行って効率的な計算を行っているものが多い。一般的に、マルチグリッド法の特徴は、単一の格子で分割するのではなく、低レベルでは細かい格子、高いレベルでは粗い格子を用いるもので、具体的にはレベルＬの格子系では、（体積ではなく辺に関して）レベル１の２のべき乗倍（２^L-1倍）のサイズの格子を用いる（図２、図３参照）。一番細かい格子は、例えば、０．１ｎｍ刻みのものである。

ポテンシャルの階層的な分離(splitting)：
Skeelらのアイデアは、いったん平滑化して分離された長距離型カーネルを、さらに、逐次的に分離するというアイデアであり、以下、具体的に説明する。図２、図３に示されている階層的な格子レベルＬに応じて、遮蔽半径を２のべき乗倍、
ｔＲ＝２^LＲ、
すなわち、Ｒ，２Ｒ，４Ｒ，８Ｒと大きくとった長距離型カーネルを導入し、次式に従って短距離型のカーネルＧ^S _tRを逐次分離していく。すなわち、
Ｇ（ｒ）＝Ｇ_R（ｒ）＋Ｇ^S _R（ｒ） （レベル０），
Ｇ_R（ｒ）＝Ｇ_2R（ｒ）＋Ｇ^S _2R（ｒ） （レベル１），
Ｇ_2R（ｒ）＝Ｇ_4R（ｒ）＋Ｇ^S _4R（ｒ） （レベル２），
Ｇ_Rmax/2（ｒ）＝Ｇ_Rmax（ｒ）＋Ｇ^S _Rmax（ｒ） （レベルＬmax−１） (3)
ここで、Ｒmaxは、最高レベルＬmaxに対する遮蔽半径の値である：
Ｒmax=２^LmaxＲ (4)
図１ｂは、長距離型カーネルの関数形を例示したものである。幾何学的には、この逐次分離は、図１ｂにおいて、１／ｒと横座標軸との間の領域をカーブＧ_R、Ｇ_2R、Ｇ_4Rなどを境として分割したことに対応する。短距離型カーネルＧ^S _tRは、具体的には、次式に従って計算されることになる：
Ｇ^S _R（ｒ）＝Ｇ（ｒ）−Ｇ_R（ｒ） （レベル０），
Ｇ^S _2R（ｒ）＝Ｇ_R（ｒ）−Ｇ_2R（ｒ） （レベル１），
Ｇ^S _4R（ｒ）＝Ｇ_2R（ｒ）−Ｇ_4R（ｒ） （レベル２），
Ｇ^S _Rmax/2（ｒ）＝Ｇ_Rmax/4（ｒ）−Ｇ_Rmax/2（ｒ） （レベル Ｌmax−１） (5)
以下では便宜上、次のような記法を採用する：
Ｇ^S _Rmax（ｒ）＝Ｇ_Rmax/2（ｒ） （レベル Ｌmax） (6)
つまり、下添字Ｒmaxには、実際には長距離(Long-range)カーネルが対応しているのだが、分離(Split)されたカーネルの意味で上添字Ｓを付加した訳である。

その結果、元のポテンシャルＧ（ｒ）は、各レベルでの短距離型のカーネルの和に分離して表される：
Ｇ（ｒ）＝Ｇ^S _R（ｒ）＋Ｇ^S _2R（ｒ）＋Ｇ^S _4R（ｒ）＋…
＋Ｇ^S _Rmax/2（ｒ）＋Ｇ^S _Rmax（ｒ） (7)
図１ｃは、このようにして得られた短距離型カーネルの関数形の例を示している。

孤立系に対する多重格子の例：
図２は、孤立系に対する多重格子の例と、それぞれのレベルの格子系に対応する短距離型カーネルを示したものである。これに応じて、（クーロンの）全相互作用エネルギーも分離される：
Ｕ＝Ｕ^S _R＋Ｕ^S _2R＋Ｕ^S _4R＋…＋Ｕ^S _Rmax/2＋Ｕ^S _Rmax (8)
ここで、
Ｕ^S _tR＝（１／２）Σ_i Σ_j ｑ_iｑ_j Ｇ^S _tR（ｒ_ij） （粒子系） (9)
マルチグリッド法の思想は、このように分離された相互作用については、カーネルの滑らかさの度合いに応じて電荷を当該レベルの格子系に粗視化して扱うのが効率的である、というものである。すなわち、レベル１以上に対しては、元の粒子の電荷を用いるのではなく、代わりに適当な方法でレベルＬの格子点に電荷を割り振っておき、格子系に対するエネルギーを次式で計算する：
Ｕ^S _tR＝（１／２）Σ_k Σ_m ｑ'_m ｑ'_k Ｇ^S _tR（ｎ_k−ｎ_m） （格子系） (10)
ここで、ｑ'_kとｑ'_mは格子点の電荷である。ｎ_k，ｎ_mは格子点の位置ベクトルである。なお、この式以降では、特に断らない限り、ｉ，ｊは粒子の番号を示し、ｋ，ｌ，ｍは格子点の番号を表すものとする。なお、レベル０における相互作用エネルギーＵ^S _Rは直接、粒子系で計算する。

周期系に対する多重格子の例：
図３は、Skeelらによる既存アルゴリズム全体を周期系に対する格子系を用いて説明した図である。図３の左側の上向きの経路（図において「(1)」で表示）は、格子点の電荷（図において黒丸で表示）を上位の粗いレベルの格子点にアサインするステップを示し、右側の下向きの経路（図において「(3)」で表示）は、格子点でのポテンシャルや電場、力（図において黒四角で表示）を下位の細かいレベルに補間内挿していくステップを示している。電荷を表す黒丸の大きさは、対応する格子点での電荷の大きさに対応し、ポテンシャルや力を示す黒四角の大きさは、対応する格子点でのポテンシャルや力の大きさに対応する。さらに図３における水平方向の経路（図において「(2)」、「(4)」で表示）は、当該レベルの相互作用カーネルを用いて、格子間相互作用を計算するステップを示している。以下の説明において、電荷を粗いレベルの格子へアサインする工程を上向き経路(Upward path)と呼び、同一レベルの格子点間に働く相互作用を計算し、格子点に作用する静電ポテンシャルを計算する工程を水平経路(Horizontal path)と呼び、粗いレベルの格子で得られた静電ポテンシャルを、逐次、細かいレベルの格子へ補間していく工程を下向き経路(Downward path)と呼ぶ。実際には、下向き経路は、ひとつ上位の粗いレベル由来の静電ポテンシャルを当該レベルの格子に補間し、そのレベルでの水平経路つまり当該レベル由来の静電ポテンシャルとの和をとる工程となっている。

分子シミュレーションに関する特許文献の例：
ＷＯ９７／４６９４９号公報（特許文献１）には、分子動力学計算を行う装置を有する分子モデル化システムが開示されている。特開２０００−２６８０６４号公報（特許文献２）には、分子シミュレーションを行う際の粒子集団の初期配置を生成する方法が開示されている。特開２００２−８０４０６号公報（特許文献３）及び特開２００３−１２５６６号公報（特許文献４）には、ポリマーやロイコ染料などの長鎖型分子の高次構造を予測するために非結合相互作用力として１２−６レナード−ジョーンズ・ポテンシャルを用いて分子シミュレーションを行うことが開示されている。特開平８−６３４５４号公報（特許文献５）には、分子シミュレーションに用いられる粒子シミュレーション方法であって、非周期境界条件と周期境界条件の両方の場合の粒子系について少ない計算量で精度よく粒子位置におけるポテンシャル値を計算できる方法として、２体ポテンシャル関数ｇの半値幅に比例する格子点を作成し、２体ポテンシャル関数ｇを２個の関数αとβとの合成積に分解し、粒子に割当てられたスカラー値を関数βにより格子点に割当て、その格子点上のスカラー値に関数αの重みを乗じる方法が開示されている。ＷＯ００／４５３３４号公報（特許文献６）には、分子動力学方法によってタンパク質の３次元モデル構造を計算することが開示されている。
ＷＯ９７／４６９４９ 特開２０００−２６８０６４ 特開２００２−８０４０６ 特開２００３−１２５６６ 特開平８−６３４５４ ＷＯ００／４５３３４ U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, "A smooth particle mesh Ewald method," J. Chem. Phys., 103(1995), 8577-8593. L. Greengard, and V. Rokhlin, "A fast algorithm for particle simulation," J. Comput. Phys., 73(1987), 325-348. R. D. Skeel, I. Tezcan, and D. J. Hardy, "Multiple grid methods for classical molecular dynamics," J. Comput. Chem. 23(2002), 673-684. W. Hockney, J. W. Eastwood, "Computer Simulation Using Particles," McGraw-Hill, New York, 1981 P. Allen and D. J. Tildesley, "Computer Simulation of Liquids," Oxford Univ. Press, 1987 E. L.Pollock and J. Glosli, "Comments on P3M, FMM, and the Ewald method for large periodic Coulombic systems," Comput. Phys. Comm., 95(1996), 93-110

しかしながら、上述したSkeelらのマルチグリッド法には、計算は速いが精度は低い、という問題点がある。特に、エネルギーの精度が悪い、という問題点を有する。そこで本発明は、以下の３つ課題を解決することを目的とする。

第１の課題は、Skeelらのマルチグリッド法の精度の低さの改良に関する。電荷ｑ_iを持つ原子ｉの近傍に格子点ｋがあってその格子点ｋの位置ベクトルをｎ_kで表すとき、Skeelらは、電荷に基底関数Ｗ（ｒ_i−ｎ_k）なる重みをつけて割り振ることにより、格子点ｋに対して電荷ｑ'_kをアサインした。ｒ_iは原子ｉの位置ベクトルである。電荷ｑ'_kのアサインを式で表わすと次式になる：
ｑ'_k＝Σ_i ｑ_i Ｗ（ｒ_i−ｎ_k） (11)
ここで、重みづけに使われた関数Ｗ（Skeelらの原著論文ではφ_kと記載）は、局所的な台(Local support)をもつ基底関数である。つまり、局所的な有限区間においてのみ関数の値が非ゼロである。この基底関数として数種の例がSkeelらの論文に記載されており、その典型的な例は区分的多項式関数の性質を有する３次のＢスプライン関数であり特定の格子点の近傍においてのみ非ゼロ値をとっている。しかし、Skeelら自身の数値実験結果の報告によれば、数種の基底関数に対して得られた計算精度は、現在必要と考えられている精度に比べて２〜３桁精度不足であり、低精度に留まっている。

分子シミュレーションの目的にもよるが、精度として、力の場合、最悪でも１×１０^-3程度の相対精度が必要であり、自由エネルギーとの比較を行う場合には、最悪でも１×１０^-4程度の相対精度が必要であるとされている。シミュレーション結果を安心して使用するためには、１×１０^-5〜１×１０^-6程度に相対精度が目安とされている。これに対し、Skeelらの結果では、最良のものの相対精度は６×１０^-4であり、かろうじて最低条件を満たしているものの、最良の結果を出すためには計算時間が多くかかっているのが現状である。Skeelらの方法は、１×１０^-3程度までの相対精度での計算ならば高速であるが、それ以上の精度が要求される場合には、ＰＭＥ法などよりも計算が遅くなるとの報告もある。言い換えれば、高い精度であっても高速で計算できるマルチグリッド法の開発が必要であると言える。

第２の課題は、孤立系に対する高速な計算方法が欠落しているという、現状を改良することである。

第３の課題は、ＰＭＥ法やＦＭＭ法はＭＤ専用ボードでは高速化が難しいことを踏まえ、専用ボードのハードウェアに適したマルチグリッド法を開発することである。

本発明は、以下の２つの手法から成り立っている。

（１）既に述べたように、Skeelらの方法ではシミュレーションで必要とされる十分な精度を達成することが不可能であった。これは、Skeelらの方法における基底関数が、１階または２階導関数にまでのみ連続であって、それ以上の高階導関数では不連続であるので、誤差が格子幅の約２乗に比例して発生しており、近似精度が悪いためであると考えられる。

本発明の第１の手法は、基底関数として少なくとも２階以上の導関数にまで連続なものを採用し、当該レベルにおいて格子点のペア（または粒子のペア）のいずれの格子点（または粒子）も、１段階上位であってより粗いレベルの格子の格子点に重なっている場合には、ひとつ上位レベルでのペアのエネルギーが必ず当該レベルの正しいエネルギー値と一致することを特徴とする、上位レベルの格子点に対する電荷決定方法を用いることをキーアイデアとする。

この手法において従来の計算ステップに追加するポイントは２箇所ある。

追加ステップのひとつは、格子点の電荷に関するものであり、いったんSkeelらの方法でひとつ上位のレベルの格子点に電荷をアサインした後、さらに広範囲の格子点に再アサインする計算ステップを追加する。図４は、これを１次元系で説明したものであり、図示される円は格子点での電荷を示し、その円の大きさは電荷の大きさを示している。（ａ）に示すように、ひとつ下位レベルの格子点（当該レベル１では原子）の座標がｘ、当該レベルの格子サイズがｈ、ｘに最も近い格子点がｋであるとき、Skeelらの方法ではｘの近傍の格子点（座標ｋｈ）のみに電荷がアサインされるが（図４の（ｂ））、本発明に基づく再アサインによって、広い範囲の格子点に電荷を再アサインされる（図４の（ｃ））。

追加ステップの２つ目は、ポテンシャルと力の内挿補間計算に関するものである。このステップは、ポテンシャル値をひとつ下位のレベルの格子点に内挿補間する前に、前記の再アサインに基づくエネルギー表式にコンシステントな方法を用いて、当該レベルの広い範囲の格子点にポテンシャルを配分するステップである。

これらの２つのステップの追加によって、本発明によれば、高い計算精度を維持することができる。

（２）第２の手法は、第１の手法を踏まえたものであり、それぞれのレベルの格子系において孤立系に対してマルチグリッド法を適用する場合の高速化に関するものである。この手法では、孤立分子系を周期ボックス（つまり、周期性の繰返し単位）に入れて、空間に周期的に配置して周期的境界条件を適用するものであり、図５はこの手法を１次元系において説明したもので、図において、丸印は、電荷を有する原子または電荷を有する格子点を示している。この手法は、異なる周期ボックスに存在する格子点が当該カーネルの影響到達距離（遮蔽距離）以上に隔離されるように、付加的な緩衝領域を孤立分子系を包含する最小のボックスの周囲に追加して周期ボックスを作るものである。このように周期性を導入することにより、３Ｄ−ＦＦＴ法を用いた計算が可能になる。付加的な緩衝領域内には、粒子は入っていない。ここで緩衝領域の辺の大きさは、影響到達距離と格子サイズから決定される。

第２の手法の利点を鮮明にするために、単一レベルのみを用いる場合と比較する。単一レベルでは、周期ボックスを互いに隔離するには、緩衝領域の辺を、孤立系を包含する最小のボックスの辺の長さ以上に設定する必要があるので、格子点数がマルチグリッド法での倍程度必要であり、それに伴って計算時間も所要メモリも大きくなる。一方、マルチグリッド法では、細かいレベルでは遮蔽半径したがって緩衝領域が小さいので、格子点数と計算時間の増加は小さく抑えられ、粗いレベルでは緩衝領域は大きいが格子点数自体が小さいので計算時間は少なくて済み、結局、計算時間が短くなる。

なお、上述した第１の手法を対象とすることなくSkeelらの方法にこの第２の手法を適用することも不可能ではないが、その場合には、十分な精度を確保できない可能性がある。第２の手法は、第１の手法と組み合わせて用いることが好ましい。

本発明の分子シミュレーション方法は、電荷が付与されている粒子を有する系における非結合相互作用をマルチグリッド法で計算する段階を有する分子シミュレーション方法において、多重格子系構築手段が、入力された原子座標に応じて多重格子系を構築して前記多重格子系をメモリに格納する段階と、対象とするレベルより１段階上位であってより粗い格子を有する上位レベルの格子点に電荷をアサインする場合に、電荷計算手段が、前記メモリにアクセスして、電荷を上位レベルでの近傍格子点にアサインし、アサインされた電荷をさらに広範囲の格子点に再アサインする段階と、対象とするレベルの格子点に作用するポテンシャル、電場または力を、対象とするレベルよりも一段階下位のより細かい格子を有する下位レベルの近傍格子点に内挿補間する際に、内挿補間に先立って、ポテンシャル計算手段が、メモリにアクセスして、対象とするレベルにおけるより広い範囲の格子点に電荷を分配する段階と、を有することを特徴とする。

このような分子シミュレーション方法において、それぞれのレベルの格子系において孤立分子系を周期系として扱うために、異なる周期ボックスに存在する格子点が対象とするレベルの格子点間相互作用カーネルの影響到達距離以上に隔離されるように、孤立分子系を包含するボックスの周囲に空の（粒子の入っていない）付加的な緩衝領域を追加することにより周期ボックスを作成する段階をさらに設け、緩衝領域のサイズは、カーネルの影響到達距離と格子サイズとから必要最小サイズを算出して決定されるようにしてもよい。周期ボックスを作成した場合には、そのような格子系に対して３次元高速フーリエ変換を適用すればよい。

本発明の分子シミュレーション装置は、電荷が付与されている粒子を有する系における非結合相互作用をマルチグリッド法で計算する段階を有する分子シミュレーション装置であって、メモリと、入力された原子座標に応じて多重格子系を構築し、メモリに格納する多重格子系構築手段と、メモリにアクセスし、対象とするレベルより１段階上位であってより粗い格子を有する上位レベルの格子点に電荷をアサインする電荷計算手段と、メモリにアクセスし、対象とするレベルの格子点に作用するポテンシャル、電場または力を、対象とするレベルよりも一段階下位でより細かい格子を有する下位レベルの近傍格子点に内挿補間するポテンシャル計算手段と、を備え、電荷計算手段は、電荷を上位レベルでの近傍格子点にアサインしてその結果をメモリに格納する手段と、メモリを参照してアサインされた電荷をさらに広範囲の格子点に再アサインしてその結果をメモリに格納する手段とを有し、ポテンシャル計算手段は、内挿補間に先立って、対象とするレベルにおけるより広い範囲の格子点に電荷を分配する。

この分子シミュレーション装置は、さらに、それぞれのレベルの格子系において孤立分子系を周期系として扱うために、異なる周期ボックスに存在する格子点が対象とするレベルの格子点間相互作用カーネルの影響到達距離以上に隔離されるように、メモリに保持された多重格子系において孤立分子系を包含するボックスの周囲に空の（粒子が入っていない）付加的な緩衝領域を追加することにより周期ボックスを作成する緩衝領域設定手段を備えていてもよい。この場合、緩衝領域のサイズは、カーネルの影響到達距離と格子サイズとから必要最小サイズを算出して決定される。

計算機が３Ｄ−ＦＦＴ計算を処理できる能力は、そのアーキテクチャに大きく依存しており、差が大きい。このような３Ｄ−ＦＦＴの処理能力の度合いに応じて、本発明によれば、それぞれ、以下に記載の効果が得られる。

ＣＰＵの台数が多い場合（現時点でのＰＣクラスタの性能で考えると約３２台以上）には、従来の方法では３Ｄ−ＦＦＴ演算が並列処理化（パラレル化）のボトルネックとなっているが、本発明によれば、３Ｄ−ＦＦＴを用いない実空間型（非ＦＦＴ型）のマルチグリッド法によって、許容精度範囲で高いパラレル化を達成することができる。

一方、ＣＰＵの台数が少ない場合（現時点でのＰＣクラスタの性能で考えて約３２台以下）には、孤立系に関しては従来は十分な高速化はできなかったが、本発明を適用することにより、３Ｄ−ＦＦＴ法が使用できることになり、許容精度範囲内で従来法よりも高速計算することが可能になる。なお、周期系でＣＰＵ台数が少ない場合に関しては、本発明のマルチグリッド法にはＰＭＥ法を超える明確な利点は存在しないので、ＰＭＥ法の使用で十分であると考えられる。

さらに本発明によれば、ＭＤ専用ボードを併用する場合に、高精度かつ高速計算が可能になる。既存法であるＰＭＥ法やＦＭＭ法は、近距離では粒子−粒子の相互作用を計算し、長距離では粒子−メッシュの相互作用を計算するものである。現在実用化されているＭＤ専用ボードは粒子―粒子の相互作用のみを高速計算するものである。粒子−メッシュ相互作用を計算するＭＤボードは、その相互作用が大変に複雑な相互作用であるために、作製自体が困難である。したがって従来は、粒子−粒子の相互作用を専用ボードでいかに高速に処理したとしても、粒子−メッシュの相互作用を汎用機で処理する部分が残ってしまい、全体としての性能を引き上げるには限界があった。本発明を用いた場合には、格子点間の短距離型カーネルを通じての相互作用は粒子間のクーロン相互作用と形式的には同じ形であるので、現在市販されている専用ボードでその短距離型カーネルを通じた相互作用が計算することが可能になり、汎用機の負担を大幅に減少させることができる。したがって、実用上十分な精度を維持しつつ計算の高速化を達成することが可能になる。

次に、本発明の好ましい実施の形態について、図面を参照して説明する。

第１の実施形態：孤立系の計算：
まず、本発明の第１の実施形態として、孤立系に本発明を適用した例を説明する。図６、図７、図９は第１の実施形態における処理の流れを示しており、特に、図６は、全体の処理の流れを示している。ここではタンパク質などの生体高分子の分子シミュレーションを対象として本発明の方法を適用した場合を説明する。図７は図６におけるステップ１０８を詳細に示したもので、その終端はステップ１０９につながる。また、図９はステップ１０９を詳細に示したもので、その終端はステップ１１０につながる。

多重格子系の構築のための手順：
まず、ステップ１０１において、生体高分子の原子座標を入力する。具体的には、対象とする生体高分子を構成する各原子の座標やその原子のタイプを記載したデータを入力する。たとえば、ＰＤＢ(Protein Data Bank)形式のデータを入力する。

次に、ステップ１０２において、原子の系に対してＡｍｂｅｒやＣｈａｒｍｍなどの分子計算プログラム（ＭＤコード）を用い、アミノ酸残基の電荷の偏在に応じて、各原子に部分電荷を割り振る。

次に、ステップ１０３において多重格子系のパラメータを設定する。具体的には、原子系を覆う格子の外枠を、原子のｘ，ｙ，ｚ座標の最大値や最小値に基づいて、若干のゆとりを含めて設定する。さらに、格子の最も粗いレベルＬmaxの値や、最も細かいレベル０における格子サイズ（つまりメッシュ幅）ｈを設定する。

次に、ステップ１０４において、各レベルの格子系において、ひとつ下位レベルの細かい格子系に基づいて格子系を構築し、これによって多重格子系を構築する。ここでは、基底関数の例として、Ｂスプライン関数を用いる。なお、レベル０は粒子系（原子系）に対応している。

よく知られているように、Ｂスプライン関数を用いて与えられた関数の補間を行った場合、区間の端点において補間精度が劣化することがある。そのため、ひとつ上位のレベルの格子系を構成するときは、Δｍ個の余分な格子を区間の両端に追加してもよい。代表的な例では、使用するＢスプラインの次数をＰとするとき、Δｍ＝（Ｐ／２−１）であり、追加しない場合はΔｍ＝０である。格子の範囲は次式で決定すればよい：
（レベルＬの下端点の座標）＝（レベルＬ−１の下端点の座標）−Δｍ（レベルＬの格子サイズ）；
（レベルＬの格子点数）＝切り上げ［（レベルＬ−１の格子点数）／２］＋２Δｍ；
（レベルＬの格子サイズ）＝（レベルＬ−１の格子サイズの２倍）＝２^L-1ｈ； (12)
ここで、格子点を（ｎ_x，ｎ_y，ｎ_z）の３整数の組で表わし、ｎ_x＝０，１，２，…，Ｍ_x；ｎ_y＝０，１，２，…，Ｍ_y；ｎ_z＝０，１，２，…，Ｍ_zの値をとっている。上記の式で格子点数とは、Ｍ_x，Ｍ_y，Ｍ_zのいずれかを指している（Ｍ_x＋１ではない）。軸方向で格子のサイズが異なるときは、適宜読み替えればよい。具体的な実施例においては、次数Ｐの値は４以上の偶数であって、代表的な値はＰ＝４〜２０である。

なお、数万時間ステップ以上にわたる分子シミュレーションの過程で、格子系の更新は、前回更新時からの原子位置変化の最大値が０．０１ｎｍ〜０．１ｎｍ程度に設定された「しきい値」を越えた場合に行うのが好ましい。

孤立系に対して３Ｄ−ＦＦＴを適用するための手順：
次に、ステップ１０５において、孤立系に対して３Ｄ−ＦＦＴ（３次元高速フーリエ変換）を適用しようとする場合には、周囲に緩衝領域を設定する。緩衝領域を設定しようとするステップは、本発明において特徴的なものである。

孤立系は周期系ではないので、周期系においてよく用いられるＦＦＴを孤立系に対して適用することはできないようにみえる。しかしながら、以下に記載の工夫を行うことで、孤立系に対してもＦＦＴを使うことができる。

孤立系がボックスの中に含まれ、ボックスにはレベル１の格子系が張られているとする。レベル１での格子サイズをｈとすると、レベルＬの格子系に対しては、格子サイズ＝２^L-1ｈ、カットオフ距離（ポテンシャルの到達距離）＝２^LＲとなる。したがって、座標軸方向に余分にΔＭ個の空の（つまり、荷電粒子のはいってない）格子を緩衝領域として追加して、２^L-1ｈΔＭ＞２^LＲ、つまり、
ΔＭ ＞２Ｒ／ｈ (13)
が成立するようなΔＭの最小値を選択してボックスを拡大しておけば、レベルＬの格子系では常に隣のボックスが見えないことになるので、周期系として扱うことができ、ＦＦＴ法を適用できることになる。

図５は緩衝領域の追加を１次元系、レベル１に関して説明したものである。短距離型のカーネルの到達距離パラメ−タＲ＝１．４、格子サイズｈ＝１．０の場合、分子系は総数Ｍ＝５の格子で覆っていたが、荷電粒子の入っていない空の格子をΔＭ＝３個、周囲に追加すると、周期的ボックスに対応する格子総数はＭ＝８に拡大される。レベル１の短距離型カーネルの到達距離は高々２Ｒ＝２．８であるので、隣のボックスとの相互作用は生じなくなる。もし、軸方向によって格子サイズが異なる場合には、軸方向ごとに異なるΔMを使用すればよい。

格子系がＭ×Ｍ×Ｍ個の格子点からなる場合には、その相互作用ペアの個数はＭ⁶個のオーダーであるので、計算量はＭ⁶のオーダーとなるが、３Ｄ−ＦＦＴ法を用いると、この計算量を、Ｍ³ ｌｏｇ Ｍのオーダーの計算量にまで減少させることができる。したがってこの方法は、３Ｄ−ＦＦＴが高速で作動できる構成のコンピュータにおいて特に有効になる。

なお、従来も孤立分子を空間内に周期的に配置することによって周期境界条件が達成されるようにすることが提案されていたが、その場合には、周期構造における分子と分子との間の間隔をその分子の大きさよりも大きくする必要があり、全体としての空間が大きくなってＦＦＴの計算点数も莫大になる、という問題が生じでいた。本実施形態の場合は、カーネルの到達距離は分子の大きさよりもかなり小さいので、従来提案されていた方法に比べ、周期構造の中で分子を十分に密に配置することが可能となり、ＦＦＴの計算点数が大幅に減って、短い時間内で計算を完了されることが可能になっている。

カーネルの分離：
次に、ステップ１０６において、各レベルにおいて、カーネルを長距離型と短距離型カーネルに分離する。

関数値表の作成：
カーネル関数値の計算に計算時間がかかるので、ステップＳ１０７において、あらかじめ短距離型カーネルの関数表を作成しておくとよい。短距離型カーネルの関数値表を作成するステップは、本発明において特徴的なものである。この実施形態では、関数値表の作成に際し、以下の方法で決定される関数を採用した。

半径Ｒの球の内部に次式に従って電荷密度ρ（ｒ）を分布させる：
ρ（ｒ）＝ｃｏｎｓｔ［（１−ｒ／Ｒ）（１＋ｒ／Ｒ）］^d-2 (14)
ただし、ｃｏｎｓｔは、規格化定数であって球内の電荷の総量が１になるように決定する。電荷密度は表面に近づくにつれて薄くなって滑らかにゼロになる。この場合の静電ポテンシャルΦ（ｒ）は、ポアソン方程式
∇²Φ（ｒ）＝−４πρ（ｒ） (15)
の解として解析的に得ることができる。原点に鋭く分布していた点電荷を、原点を中心とする滑らかな電荷分布に置き換えたので、式(15)の解として得られる静電ポテンシャルは、原点から距離Ｒ以上離れた遠方では１／ｒに等しいが、原点に近づいても発散せずに有限値に収束するので、この解は長距離型カーネルＧ_R（ｒ）として相応しい性質を有する。関数の滑らかさに関して重要なポイントは、電荷分布がｄ−２回連続微分可能なので、Ｇ_R（ｒ）はｒ＝Ｒでｄ回連続微分可能になることである。スプライン関数による補間誤差は、スプラインの次数をＰとするとＧ^S _R（ｒ）のＰ階微係数によって上限が定まるので、スプラインの次数Ｐと、このｄとを一貫させることは、最小のコストで精度を維持するために重要である。

そのほかにも、以下に列挙したような様々な関数を用いて分離することが可能である：
(1)Skeelらの文献に記載の関数例（非特許文献３のFig. 3；文献中のａが本明細書の式(2)のＲに相当），
(2) Hockneyらの文献で定義された関数（非特許文献４のeq(8-75)；文献中のａが本明細書の式(2)のＲに相当），
(3)誤差関数型の関数（非特許文献５に記載のEwald和）：
Ｇ_R（ｒ）＝［１−ｅｒｆｃ（αｒ）］／ｒ ，
Ｇ_R（０）＝２α／π^1/2，
Ｇ^S _R（ｒ）＝ｅｒｆｃ（αｒ）／ｒ (16)
ここで、ｅｒｆｃは、相補的誤差関数(complementary error function)のことである。１／α＝３Ｒと対応づけると、ｒ＞ＲでＧ^S _R（ｒ）は（完全ではないが）実質的にゼロになる。

粒子に作用するポテンシャルと力の計算：
短距離カーネルの関数値表が作成されたら、次に、この関数値表を用いて、粒子に働くポテンシャルと力の計算を開始する。ここでは、当該レベルの原子ｊが、短距離型カーネルを通じて他の原子から受けるポテンシャルを次式にしたがって計算するものとする。

φ^S（ｒ_i）＝Σ_j ｑ_j Ｇ^S _R（ｒ_i−ｒ_j） (17)
ここで、ｑ_jは、原子ｊの電荷である。これは、ごく普通の粒子−粒子相互作用である
まず、ステップ１０８において、レベル０（粒子系）に対して、短距離型カーネルに基づき、粒子に働くポテンシャルと力を計算する。これは上述の図３において(4)で示した水平経路の計算に相当する。

上向き経路での計算（格子レベルＬを１からＬmaxまで逐次増加）：
〈上向き経路の計算の開始〉
次に、ステップ１０９において、格子上の電荷を計算するために、上述の図３において説明した上向き経路(Upward path)の計算を開始する。このパス（経路）では、格子点の電荷を逐次決定してゆく。

図７は、上向き経路の計算をフローチャートとして示したものである。ここに示すように上向き経路の計算は、要約すれば、上述したステップ１０８においてレベル０に対して短距離型カーネルに基づき粒子に作用するポテンシャルと力とを計算したのち、レベルＬが１からＬmaxまで１ずつ増加するループを開始し（ステップ１２１）、このループにおいて、ひとつ下位のレベルの格子点に基づいて、当該レベルの格子点に電荷をアサインし（ステップ１２２）、いったんアサインされた電荷に基づいて、広範囲の格子点に電荷を再アサインし（ステップ１２３）、フーリエ空間（ＦＦＴ）で計算かもしくは実空間（非ＦＦＴ）で計算かを選択し（ステップ１２４）、この選択にしたがって各格子点におけるポテンシャルを、当該レベルＬでの短距離型カーネルに基づいて計算し（ステップ１２５）、その後、ステップ１２１で開始したループの終端に至る（ステップ１２６）というものである。ステップ１２３の再アサインを行うステップは、本発明における特徴的なステップである。ステップ１２２は、上述した図３において(1)で示した上向き経路の計算に対応し、ステップ１２５は、上述した図３において(4)で示した水平経路の計算に対応する。

以下、上向き経路の計算について、当該レベルＬ＝１の格子系の場合に対して説明する。

この場合、下位レベル０が粒子系になり、ｔ＝２^Lは２に等しい。この格子点の電荷を、粒子が格子点に一致する場合には正しいエネルギーと一致するように決定すれば、粒子系のエネルギーを格子系で精度良く近似できることになる。

まず数学的背景を説明する。カーネルＧ^S _tR（ｒ_i−ｒ_j）をｒ_iとｒ_jの関数とみなして基底関数Ｗを用いて２重展開し、ｒ_iとｒ_jがともに格子点に一致する場合には必ず、元の関数値（つまり２粒子系に対する正しいエネルギー値）に一致するという条件を課して展開係数を決定すると、以下の展開近似式（補間式）が得られる。

Ｇ^S _tR（ｒ_i−ｒ_j）＝Σ_k Σ_m Ｗ（ｒ_i−ｎ_k）Ｗ（ｒ_j−ｎ_m）×Σ_k' Σ_m'（Ａ^-1）_kk'（Ａ^-1)_mm'Ｇ^S _tR（ｎ_k'−ｎ_m'） (18)
Ｂスプラインを基底関数Ｗに用いる場合の展開誤差は、Ｇ^S _tRのＰ階導関数の絶対値の上限とｈ^Pとの積のオーダーであるので、補間式は精度の高いものになる。この展開近似式を次式で与えられるエネルギー表式：
Ｕ^S _tR＝（１／２）Σ_i Σ_j ｑ_i ｑ_j Ｇ^S _tR（ｒ_i−ｒ_j） （粒子系） (19)
に代入して得られたのが、格子系に対するエネルギー表式である：
Ｕ^S _tR＝（１／２）Σ_k,l ｑ'_k ｑ'_lＧ^S _tR（ｎ_k−ｎ_l） （格子系） (20)
この式に含まれている格子点上の電荷ｑ'_kは、以下の手続きで計算される。

〈電荷のアサイン（assign；割り振り）ステップ〉
荷電粒子の近傍の格子点に重みＷ（ｒ_j−ｎ_k）をつけて電荷をアサインする（図４の（ｂ）を参照）：
ｑ"_k＝Σ_j ｑ_j Ｗ（ｒ_j− ｎ_k） (21)
この式で、ｎ_kは格子点ｋの位置ベクトルであり、また、粒子ｊに関する和は、Ｗ（ｒ_j−ｎ_k）がゼロと異なる値をとる粒子についてとればよい。それらの粒子ｊは格子点ｋの近傍に位置している。なお、Skeelらによって提案されている従来のマルチグリッド法は、このアサインのステップのみで次の電荷の再アサインステップが存在しないため、低精度である。

ここに現れているＷ関数は、通常のＢスプライン関数の積であり、次式のように、ｘ，ｙ，ｚ成分にデカップル(decouple)した形で表される：
Ｗ（ｒ_j−ｎ_k）＝Ｂ（（ｒ_jx−ｎ_kx）／ｈ）Ｂ（（ｒ_jy−ｎ_ky）／ｈ）×Ｂ（（ｒ_jz−ｎ_kz）／ｈ） (22)
Ｂスプライン関数の定義・性質については多く教科書で説明されているが、その中で、以下の説明で重要な点は、Ｐ次のＢスプライン関数Ｂ（ｘ）は、区間（−Ｐ／２，＋Ｐ／２）において区分的Ｐ−１次多項式であり、区間外ではゼロをとるという性質である。従って、粒子ｊに関する和は、たとえば、ｘ成分に関しては、格子点ｋの近傍に位置し｜ｒ_jx−ｎ_kx｜／ｈがＰ／２以下の粒子ｊについてとれば十分である。

〈電荷の再アサインのステップ〉
次式にしたがって、広い範囲の格子点に電荷を再アサインする（図４の（ｃ））：
ｑ'_l＝ Σ_k ｑ"_k（Ａ^-1）_kl (23)
ここで、Ａは、次式で定義される行列であり、いわば、重み行列Ｗの逆行列である：
Ａ_kl＝Ｗ（ｎ_k−ｎ_l） (24)
なお、上記の式から分かるように、格子点の電荷値はカーネルの具体的な関数形には依存していない。

関数Ｗはｘ，ｙ，ｚ成分にデカップルされているので、行列ＡならびにＡの逆行列もｘ，ｙ，ｚ成分にデカップルされて、Ａ^-1＝Ａｘ^-1 Ａｙ^-1 Ａｚ^-1の形に表される。したがって、格子点ｋを３次元で顕に（ｋｘ，ｋｙ，ｋｚ）と書くと、上記の式のｋについての和は、この３成分を順番に実行すればよい：
ｑ'_lx,ly,lz＝Σ_kz｛Σ_ky［Σ_kx ｑ"_kx,ky,kz（Ａｘ^-1）_kx,lx］（Ａｙ^-1）_ky,ly｝（Ａｚ^-1）_kz,lz (25)
Ｗ（ｎ_k−ｎ_l）がゼロでない値をとるのは２つの格子点が空間的に近接する場合のみであり、その近接の度合いはスプラインの次数Ｐに応じて定まっているので、行列Ａもそのｘ，ｙ，ｚ成分も実際にはバンド行列となり、バンド行列連立方程式に対する標準的な解法で解くことができる。上記の行列Ａやその逆行列の更新は、格子系が再構築された場合に行えば十分である。

図８ａは、スプライン次数Ｐ＝６に対して、行列Ａｘ，Ａｙ，Ａｚの例を示したものである。追加された余分な格子点がｋｘ＝−２，−１である。なお、格子の終端でのＡの成分は、始端での値と点対称の形になる。

上記の説明は、当該レベルＬ＝１に対するものであった。高位のレベルに関しては上記の説明を、粒子は「ひとつ下位の細かいレベルの格子点」、格子点は「当該レベルの格子点」で読み替えればよい。

〈高階の導関数値を０とおく補間方法を採用する場合の変更点〉
端点での高階微係数をゼロとする境界条件を採用する場合には、行列Ａの表式は端点で変更されることになる。取り扱いは複雑になるが、若干ではあるが、精度向上や、行列サイズ減少に伴う計算処理速度向上という利点がある。

変更点を以下に説明する。格子系構築の際に、Δｍ＝Ｐ／２−１個の余分な格子点を区間の両端に追加したものとする。ここでは、両端にΔｍ個の格子点を追加する以前に存在している格子点をｋとｍで表し、その中でも特に両端に位置していた格子点をｅ（endsの意）で表し、また、両端に追加されたΔｍ個の格子点に関してはａ（additionalの意）で表す。カーネルの高階導関数の端点での補間値のうち、微分階数ｄがＰ／２次から（Ｐ−２）次までのものの値をゼロと置くという条件を付加することにする。

条件を付加した場合には、追加した格子点（つまりａ）も含めて全ての格子点の電荷を決める連立方程式の係数行列Ａ'の逆行列の計算が、式(24)に与えられるＡに代わって必要になる。

ここで、行列α，β，γ，εの成分は次式で与えられる：
α_km＝Ｗ（ｎ_k−ｎ_m），
β_ka＝Ｗ（ｎ_k−ｎ_a），
γμ_k＝Ｗ^(d)（ｎ_e−ｎ_k），
εμ_a＝Ｗ^(d)（ｎ_e−ｎ_a） (27)
ここで、μは、ｅとｄのペアを１文字で表記したものである。しかし、この行列Ａ'のバンド幅は、行列内部でＰ／２−１であるが、端ではＰ−２に広がっている。たとえば、図８ｂは行列Ａ'のｘ成分の例を示しており、ｄ＝４の行（ｅ＝始端）ではバンド幅は４に広がっており、数値計算に時間がかかる不都合がある。

しかし、両端に追加された格子点（つまりａ）の電荷については、式(27)の連立方程式の一部を解析的に解くことができるので、Ａ'からｋｍ成分以外の成分を消去することが可能であり、次式で定義される行列Ａを用いれば良いことが示される：
Ａ＝α＋βε^-1γ (28)
このＡに対しては、バンド幅は、Ｐ／２−１に収まり、好都合である。このＡを用いる場合には、追加した２Δｍ個の格子点は、定式化のために補助的に導入した仮想的なものとなり、実際に電荷を付与する必要はないし、格子点に関する和においても除外できる。図８ｃは、この行列Ａのｘ成分の例を示したものである。

〈格子点におけるポテンシャルの計算（水平経路に対応）〉
当該レベルの格子点が、短距離型カーネルを通じて他の格子点から受けるポテンシャルを、以下に記載の実空間法もしくはフーリエ(Fourier)空間法のいずれかを選択して計算する。

（１）実空間法（３Ｄ−ＦＦＴを使用せず）の場合：
次式にしたがって計算する。

φ^S（ｎ_k）＝Σ_m ｑ'_m Ｇ^S _tR（ｎ_k−ｎ_m） (29)
ここで、ｑ'_mは、格子点ｍの電荷である。カーネルの影響の到達距離は高々ｔＲであるので、ｍについての和は、格子点ｋから距離ｔＲ内に位置する格子点に対してとれば十分である（ｔ＝２^L）。

（２）フーリエ空間法（３Ｄ−ＦＦＴを使用）の場合：
よく知られているように（たとえば、非特許文献４）、格子系が周期的に並んでいる系に対するエネルギー、ポテンシャル、電場は、以下のスキームに従って、３Ｄ−ＦＦＴを用いて高速計算できる：

ここで、下添字のＦはフーリエ変換された量であることを示しており、＊は畳み込み(Convolution)積である。また、Ｇ^S _F（ｋ_x，ｋ_y，ｋ_z）は、Ｇ^S（ｎ_x，ｎ_y，ｎ_z）の（有限）フーリエ変換であって、変換の際の整数ｎ_xの変域は０からＭ_xではなく、−∞から＋∞に拡張したものである。しかし、Ｇ^S（ｎ_x，ｎ_y，ｎ_z）がゼロでない値をとるのは、ｎ＝｜ｎ｜＜２^LＲの時のみであるので、その有限範囲内で和をとれば十分である。

以上の処理によって、水平経路の計算が終了する。

下向き経路の計算（格子レベルＬをＬmax−１から０まで逐次減少）：
〈下向き経路の計算の開始〉
次に、図６のステップ１１０において、格子点でのポテンシャル、電場、力を逐次決定していくために、上述の図３において説明した下向き経路(Downward path)の計算を開始する。

前述の電荷の再アサインメントで得られたエネルギー表式は、厳密値とは異なる近似的なものであるが、静電ポテンシャルや力は、一貫して、この近似的なエネルギー表式を粒子の座標で微分して得られた計算式を用いるものとする。この意味で、電荷の再アサインとコンシステントな静電ポテンシャルや力と呼ぶことにする。以下では、それらの計算式に基づいて、ひとつ上位の粗い格子で得られたポテンシャル値から、当該レベルの格子点のポテンシャル値を計算してゆく。

図９は、下向き経路の計算をフローチャートとして示したものである。ここに示すように下向き経路の計算は、要約すれば、レベルＬがＬmaxから０まで１ずつ減少するループを開始し（ステップ１３１）、このループにおいて、ひとつ上位のレベル（Ｌ＋１）の格子点でのポテンシャル値を、電荷の再アサイメントとコンシステントな方法でレベルＬ＋１の格子に分配し（ステップ１３２）、分配されたポテンシャル値に基づいて、レベルＬの格子点のポテンシャル値を補間し（ステップ１３３）、格子点においてレベルＬに対して計算されたポテンシャル値と、上位レベルの格子からの内挿補間を経て得られたポテンシャル値との和として、格子点におけるポテンシャルを計算し（ステップ１３４）、その後、ステップ１３１で開始したループの終端に至る（ステップ１３５）というものである。ステップ１３２でのポテンシャル値を分配するステップは、本発明において特徴的なステップである。ステップ１３３は、上述した図３において(3)で示した下向き経路の計算に対応し、ステップ１３４は、上述した図３において(3)で示した下向き経路の計算値と(4)で示した水平経路の計算値との和に対応する。

以下、下向き経路の計算について、詳細に説明する。

〈ポテンシャルの分配のステップ〉
このステップは、上向き経路において電荷の再アサインを行ったことに伴って派生的に生じたものである。ここでは、当該レベルＬのひとつ上位レベル（親）の粗い格子系において、格子点のポテンシャルを他の格子点に次式に従って分配する：
φ"^S（ｎ_m）＝Σ_k （Ａ^-1）_mkφ^S（ｎ_k） (31)
なお、Skeelらの方法は、このステップを欠いており、以下の式でφ"=φと置いたものと等価である。

〈ポテンシャルの内挿補間のステップ〉
分配されたポテンシャル値φ"に基づいて、当該レベルＬの格子点ｉ（または粒子ｉ）に働くポテンシャルφ^Sを次式のように内挿補間する：
φ^S（ｒ_i）＝Σ_m Ｗ（ｒ_i−ｎ_m）φ"^S（ｎ_m） (32)
もし、必要ならば、電場Ｅ^S、力Ｆ^Sも分配されたポテンシャル値を使って同様に内挿補間できる。ここで電場Ｅ^S、力Ｆ^Sは、いずれもベクトル量である：
Ｅ^S（ｒ_i）＝―Σ_m ∇Ｗ（ｒ_i−ｎ_m）φ"^S（ｎ_m） (33)
Ｆ^S _i＝ｑ_i Ｅ^S（ｒ_i） (34)
ここで、ｒ_iは粒子ｉの位置ベクトルであり、ｎ_mは格子点ｍの位置ベクトルである。∇は位置ベクトルに関する微分を表す。∇Ｗの計算に関しては、Ｐ次のＢスプライン関数の場合には、微分は２個のＰ−１次のＢスプライン関数を使って簡単に表わされるので、容易に計算することができる。

〈ポテンシャルの和の計算のステップ〉
上位レベルの格子から分配・補間で得られたポテンシャル値と、当該レベルに関して格子点間の短距離型カーネルを通じて得られたポテンシャル値との和を計算する。すなわち：
（当該レベル以上に由来するポテンシャル）＝（親以上のレベル由来のポテンシャル）＋（当該レベル由来のポテンシャル） (35)
を計算する。

以上の処理によって、下向き経路の計算が終了する。

結果の出力：
その後、図６のステップ１１１において、下向き経路の計算によって求められた、各粒子（原子）に作用するポテンシャルと力とを出力して、非結合相互作用の計算が終了する。

計算例：
以下、本実施形態に基づくアルゴリズムによって実際に計算を行った例を説明する。

本アルゴリズムでは、計算時間と精度は、一番細かいメッシュのサイズ（ｈ）、ポテンシャルを直接計算部とメッシュ計算部に分離する際の遮蔽距離Ｒ、スプライン次数Ｐ、および粒子位置と電荷により制御される。表１に約５万原子の系でＰ＝６に対する適用例を示す。アミロマルターゼというタンパク質とその周囲に存在する水とからなる系について非結合相互作用のシミュレーション演算を行った。シミュレーションに用いたＣＰＵは、Ｐｅｎｔｉｕｍ（登録商標）４ （動作クロック周波数＝２ＧＨｚ）である。なお同じパラメータを用いて、本発明の方法によらずに直接的な計算を行った場合には、計算に７００秒を要した。

比較のために、同じ系に対するSkeelらの計算結果（非特許文献３）を引用しておくが、かれらの精度の最も良い結果でも、力の相対誤差は、０．００２であって、本発明のものはSkeelらの結果に比べて３桁程度良好であった。

分子シミュレーション装置：
次に、上述した手順に基づいて非結合相互作用を計算する分子シミュレーション装置を説明する。図１０はこのような分子シミュレーション装置を示している。

分子シミュレーション装置は、入力部１１、多重格子系構築部１２、メモリ１３、緩衝領域設定部１４、カーネル計算部１５、電荷計算部１６及びポテンシャル計算部１７を備えている。入力部１１には、対象とする分子の原子座標や計算に必要なパラメータを入力する。多重格子系構築部１２は、上述したステップ１０２〜１０４の処理を実行することによって、原子に部分電荷を付加し、多重格子系のパラメータ設定を行い、上述した多重格子系を構築する。構築された多重格子系はメモリ１３に格納される。メモリ１３は、非結合相互作用の計算の際のワークメモリとして機能するものであって、多重格子系について、格子点ごとの電荷や電場、ポテンシャル、力などの情報も格納する。緩衝領域設定部１４は、メモリ１３に格納されている多重格子系に対し、上述したステップ１０５の処理を実行することによって、必要に応じて、その周囲に緩衝領域を設定する。カーネル計算部１５は、ステップ１０６のカーネルの分離と、ステップ１０７の関数値表の作成の処理を実行する。作成された関数値表もメモリ１３に格納される。

電荷計算部１６は、メモリ１３にアクセスして、ステップ１０８のレベル０におけるポテンシャルと力の計算と、ステップ１０９の上向き経路の計算を行う。この電荷計算部１６には、レベル０におけるポテンシャルと力を計算する初期値計算部２０と、近接格子点のペアｎ_k，ｎ_lと基底関数Ｗに対して、行列要素Ａkl＝Ｗ（ｎ_k−ｎ_l）をあらかじめ計算し、得られたバンド行列Ａの非零要素をメモリ１３に格納する第１の計算部２１と、各粒子ｑ_jの電荷に重みＷ（ｒ_j−ｎ_k）を乗じて、粒子近傍の格子点ｎ_kに割り振ることにより、各格子点上の単純アサイン電荷ｑ_k"を計算し、これをメモリ１３に格納する第２の計算部２２と、メモリ１３に格納された単純アサイン電荷ｑ_k"に対し、連立方程式Σ_mＡ_kmｑ_m'＝ｑ_k"を解くことにより、格子点に再アサインされる電荷ｑ_m'を決定する第３の計算部２３と、を備えている。第１の計算部２１は、上述の式(24)または(28)に対応する計算を行うものであり、第２の計算部２２は式(21)に対応する計算をするものであり、第３の計算部２３は、式(23)に対応する計算をするものである。

ポテンシャル計算部１７は、メモリ１３にアクセスして、ステップ１１０の下向き経路の計算を行う。

第２の実施形態：周期系での計算：
次に、本発明の第２の実施形態を説明する。ここでは、同一レベルの格子間で相互作用を計算することにより、周期境界条件を用いるがＦＦＴを使用せず、ＭＤ専用ボードでの実行に適した計算方法を説明する。

第１の実施形態の場合と同様に、対象とする分子の原子座標を入力し、各原子に部分電荷を付加し、多重格子系のパラメータを設定する。その後、多重格子系を構築する。

周期系の場合の多重格子系の構成方法
本実施形態で対象としている周期系の場合には、余分に格子を追加する必要はない。すなわち緩衝領域を導入する必要はない。しかし、採用する格子レベルの最大値をＬmaxとすると、Ｌmaxまでの全てのレベルの格子系が周期の単位ボックスとぴったりとマッチするには、（レベル０の格子点数）＝（２^Lmaxの整数倍）、 となる必要が有る。ここでも、格子点数とは、Ｍ_x，Ｍ_y，Ｍ_zのいずれかを指している。

相互作用の計算：
多重格子系を構築した後、各相互作用を計算する。従来のＰＭＥ法やＦＭＭ法では、粒子−粒子間相互作用計算部分のみを専用ハードウェアで処理しているものの、他の部分はアルゴリズムが複雑すぎるために、専用ハードウェアでの処理が実現していなかった。本実施形態によれば、格子−格子相互作用の計算にもＭＤ計算用専用ボード（専用ハードウェア）を用いることが可能になり、計算加速上、大きな効果がある。これが精度良い計算結果を与えるのは、格子点での電荷を再アサインによってうまく決定した仕組みによる。

専用ボードとは、簡単に言えば、（ｄＧ^S _tR／ｄｒ）／ｒの関数表をあらかじめ作成して格納しておき、まず、中心粒子として少数のｉ粒子の座標と電荷を演算部に読み込み、次に、他の多数の粒子（ｊ粒子）を逐次読み込んで、原子ｉに働くポテンシャルや力をパラレルに計算・蓄積していく演算装置である。

マルチグリッド法では格子レベルが複数個存在するので、専用ボードを用いる場合であっても、関数表はレベルごとに作成しておく必要がある。

各原子の感じるポテンシャルφ（ｎ_k）を計算する機能が専用ボードに備わっている場合には、上述した第１の実施形態に記載の補間法を用いて、原子に作用するポテンシャルや電場、力を直ちに計算できる。しかしながら現在入手可能な専用ボードでは、各原子の感じるポテンシャルφ（ｎ_k）を計算するための回路は省略されており、マルチグリッド法を、そのまま適用することは不可能である。その場合の解決策としては、以下の手順によって、ポテンシャル値を経由せずに電場を計算することが挙げられる：
（１）全粒子の座標と電荷（ｘｙｚｑ）のリストが座標メモリに格納されるが、このリストを２重化し、座標メモリの前半部分ではリストに本来の電荷が格納され、後半部分ではすべての電荷を１に設定する；
（２）相互作用の計算に際し、中心粒子ｉは座標メモリ後半部分から読み込み、相手となる粒子ｊは座標メモリ前半から読み込む；
（３）中心粒子ｉに作用する電場を計算する。このとき、実際に計算されるのは力であるが、ｑ_i＝１と仮定しているので、電場に等しくなる；
（４）以下に記載する微分補間型の補間法を用いて、各粒子に働く電場や力を計算する。

微分補間型の補間法
上述の第１の実施形態での補間法では、補間→微分の順に処理を行った。一方、以下で述べるのは、微分→補間の順に行う微分補間の処理に関するものである。

専用ボードを用いて、ひとつ上位の粗いレベルの各格子点ｎ_kに働く電場Ｅ^S（ｎ_k）（ベクトル量）を計算する：
E^S（ｎ_k）＝−Σ_m ｑ'_m∇Ｇ^S（ｎ_k−ｎ_m） (36)
ここで、ｑ'_mは格子点ｍの電荷を表している。

その後、汎用型のコンピュータを用いて、ひとつ上位レベルにおいて、電場の値を他の格子点に分配する。すなわち、
Ｅ"^S（ｎ_m）＝Σ_k （Ａ^-1）_mk Ｅ^S（ｎ_k） (37)
とする。そして、汎用型のコンピュータを用いて、ひとつ上位レベルで分配されて得られた電場を、当該レベルの格子点ｉ（または粒子ｉ）に内挿補間する。すなわち、
Ｅ^S（ｒ_i）＝Σ_m Ｗ（ｒ_i−ｎ_m） Ｅ"^S（ｎ_m） (38)
Ｆ^S _i＝ｑ_i Ｅ^S（ｒ_i） (39)
とする。これにより、周期系において、ＦＦＴを用いることなく、格子点あるいは粒子における電場を求めることができる。

なお、この微分補間型の補間では、上記の力の計算式に厳密に対応する「全エネルギー表式」が明確ではない（解析的な式では書き下せない）という問題点があり、ＭＤ計算中のエネルギー保存も若干悪くなるが、実用上は問題になることはない。

上述した各実施形態における非結合相互作用の計算は、一般に、それを実行するためのプログラムをコンピュータに読み込ませ、そのプログラムを実行させることによって実行される。したがって、本発明の範疇には、上述したような非結合相互作用の計算を実行するためのプログラムや、そのようなプログラムを格納した記録媒体も含まれる。

相互作用ポテンシャル（カーネル）の多重分離を説明する図であって、長距離型と短距離型への分離を説明する図である。 相互作用カーネルの多重分離と長距離型カーネルを説明する図である。 相互作用カーネルの多重分離と短距離型カーネルを説明する図である。 多重格子系と短距離型カーネルを説明する図である。 マルチグリッド法の既存アルゴリズムを説明する図である。 電荷のアサインメントと再アサインメントを１次元系において模式的に示した図である。 孤立系を周期系として扱うために孤立系の周囲に設定した緩衝領域を１次元系で模式的に示した図である。 第１の実施形態での処理を示すフローチャートである。 第１の実施形態での上向き経路の処理を示すフローチャートである。 電荷の再アサインメントに使用する行列Ａの例を示す図である。 電荷の再アサインメントに使用する行列Ａであって、格子系の端点において高階導関数をゼロとおいた場合の例を示す図である。 電荷の再アサインメントに使用する行列Ａであって、格子系の端点において高階導関数をゼロとおき、バンド幅縮小の工夫を行った場合の例を示す図である。 第１の実施形態での下向き経路の処理を示すフローチャートである。 第１の実施形態を実行するための分子シミュレーション装置のブロック図である。

符号の説明

１１ 入力部
１２ 多重格子系構築部
１３ メモリ
１４ 緩衝領域設定部
１５ カーネル計算部
１６ 電荷計算部
１７ ポテンシャル計算部
２０ 初期値計算部
２１ 第１の計算部
２２ 第２の計算部
２３ 第３の計算部

Claims (7)

1. 電荷が付与されている粒子を有する系における非結合相互作用をマルチグリッド法で計算する段階を有する分子シミュレーション方法において、
   多重格子系構築手段が、入力された原子座標に応じて多重格子系を構築して前記多重格子系をメモリに格納する段階と、
   対象とするレベルより１段階上位であってより粗い格子を有する上位レベルの格子点に電荷をアサインする場合に、電荷計算手段が、前記メモリにアクセスして、前記電荷を前記上位レベルでの近傍格子点にアサインし、前記アサインされた電荷をさらに広範囲の格子点に再アサインする段階と、
   前記対象とするレベルの格子点に作用するポテンシャル、電場または力を、対象とするレベルよりも一段階下位のより細かい格子を有する下位レベルの近傍格子点に内挿補間する際に、前記内挿補間に先立って、ポテンシャル計算手段が、前記メモリにアクセスして、前記対象とするレベルにおけるより広い範囲の格子点に前記電荷を分配する段階と、
   を有することを特徴とする分子シミュレーション方法。
2. それぞれのレベルの格子系において孤立分子系を周期系として扱うために、異なる周期ボックスに存在する格子点が前記対象とするレベルの格子点間相互作用カーネルの影響到達距離以上に隔離されるように、緩衝領域設定手段が、前記メモリにアクセスして、前記孤立分子系を包含するボックスの周囲に、前記粒子が入っていない空の付加的な緩衝領域を追加することにより周期ボックスを作成する段階をさらに有し、
   前記緩衝領域のサイズは、カーネルの影響到達距離と格子サイズとから必要最小サイズを算出して決定される、請求項１に記載の分子シミュレーション方法。
3. 前記周期ボックスが作成された格子系に対して３次元高速フーリエ変換を適用する段階を有する、請求項２に記載の分子シミュレーション方法。
4. 電荷が付与されている粒子を有する系における非結合相互作用をマルチグリッド法で計算する段階を有する分子シミュレーション装置であって、
   メモリと、
   入力された原子座標に応じて多重格子系を構築し、前記メモリに格納する多重格子系構築手段と、
   前記メモリにアクセスし、対象とするレベルより１段階上位であってより粗い格子を有する上位レベルの格子点に電荷をアサインする電荷計算手段と、
   前記メモリにアクセスし、前記対象とするレベルの格子点に作用するポテンシャル、電場または力を、対象とするレベルよりも一段階下位のより細かい格子を有する下位レベルの近傍格子点に内挿補間するポテンシャル計算手段と、
   を備え、
   前記電荷計算手段は、電荷を前記上位レベルでの近傍格子点にアサインしてその結果を前記メモリに格納する手段と、前記メモリを参照して前記アサインされた電荷をさらに広範囲の格子点に再アサインしてその結果を前記メモリに格納する手段とを有し、
   前記ポテンシャル計算手段は、前記内挿補間に先立って、前記対象とするレベルにおけるより広い範囲の格子点に前記電荷を分配する、分子シミュレーション装置。
5. それぞれのレベルの格子系において孤立分子系を周期系として扱うために、異なる周期ボックスに存在する格子点が前記対象とするレベルの格子点間相互作用カーネルの影響到達距離以上に隔離されるように、前記メモリに保持された多重格子系において前記孤立分子系を包含するボックスの周囲に前記粒子が入っていない空の付加的な緩衝領域を追加することにより周期ボックスを作成する緩衝領域設定手段をさらに備え、
   前記緩衝領域のサイズは、カーネルの影響到達距離と格子サイズとから必要最小サイズを算出して決定される、請求項４に記載の分子シミュレーション装置。
6. 分子シミュレーションの実行に際して、電荷が付与されている粒子を有する系における非結合相互作用をマルチグリッド法で計算するコンピュータを、
   入力された原子座標に応じて多重格子系を構築し、メモリに格納する多重格子系構築手段、
   前記メモリにアクセスし、対象とするレベルより１段階上位であってより粗い格子を有する上位レベルの格子点に電荷をアサインする電荷計算手段であって、前記電荷を前記上位レベルでの近傍格子点にアサインしてその結果を前記メモリに格納する手段と、前記メモリを参照して前記アサインされた電荷をさらに広範囲の格子点に再アサインしてその結果を前記メモリに格納する手段と、備える電荷計算手段、
   前記メモリにアクセスし、前記対象とするレベルの格子点に作用するポテンシャル、電場または力を、対象とするレベルよりも一段階下位のより細かい格子を有する下位レベルの近傍格子点に内挿補間するポテンシャル計算手段であって、前記内挿補間に先立って、前記対象とするレベルにおけるより広い範囲の格子点に前記電荷を分配するポテンシャル計算手段、
   として機能させるプログラム。
7. 前記コンピュータを、さらに、
   それぞれのレベルの格子系において孤立分子系を周期系として扱うために、異なる周期ボックスに存在する格子点が前記対象とするレベルの格子点間相互作用カーネルの影響到達距離以上に隔離されるように、前記メモリを参照して、前記メモリに保持された多重格子系において前記孤立分子系を包含するボックスの周囲に、前記粒子が入っていない空の付加的な緩衝領域を追加することにより周期ボックスを作成する緩衝領域設定手段であって、カーネルの影響到達距離と格子サイズとから必要最小サイズを算出して前記緩衝領域のサイズを決定する緩衝領域設定手段として機能させる、請求項６に記載のプログラム。

Images