(57)【要約】【課題】多体問題計算装置のハードウェア量を削減しつつ、計算時間の増大を防止する。【解決手段】力総和計算部２９は、基本セル内の粒子を順次計算対象粒子とし、その粒子に対して働く力の総和をイメージセル単位で求めている。ビリアル計算部３０は、系全体のビリアルを求めるための処理を、計算量のオーダーは大きいが乗算は不要な部分と、計算量のオーダーは小さいが乗算が必要になる部分とに分割した際の、計算量のオーダーが大きくなる部分についての計算を、力総和計算部２９の計算結果を利用して実行する。ソフトウェア実行部３は、ビリアル計算部３０の計算結果を利用して、計算量のオーダーが小さくなる部分についての計算を実行し、系全体のビリアルを計算する。【選択図】図３
2010072680.tif 000002
【特許請求の範囲】
【請求項１】
周期的境界条件を適用する系における基本セル内の粒子の１つを計算対象粒子として選択し、前記周期的境界条件のイメージセルを順次処理単位にして、そのイメージセル内の粒子が前記計算対象粒子に対して及ぼす力の総和を求めるという処理を、前記基本セル内の全ての粒子を計算対象粒子として選択するまで繰り返し行うと共に、前記基本セル内の粒子毎にその粒子に働く力の総和を求めて保持する力総和計算部と、
該力総和計算部においてイメージセル単位の力の総和が求められる毎に、その力の総和の各座標成分を順次選択する選択回路と、
前記力総和計算部において処理単位にされているイメージセルと前記基本セルとの位置関係に基づいて、前記選択回路で選択された力の総和の座標成分をそのまま出力するか、符号を反転して出力するか、或いは出力を０とする加減算制御回路と、
ビリアルの成分毎の加算結果が保持される保持回路と、
前記加減算制御回路の出力と、前記保持回路に保持されている加算結果の内の、前記加減算制御回路の出力に対応する成分の加算結果とを加算し、加算結果を前記保持回路に格納する加算器と、
ビリアルの成分毎に、前記基本セル内の粒子の粒子座標と前記力総和計算部に保持されている前記基本セル内の粒子に働く力の総和とを用いて、その成分についての力の総和と粒子座標との積和を求める第１の処理と、前記保持回路に保持されているその成分についての加算結果と定数とを乗算する第２の処理と、前記第１の処理で求めた積和結果と前記第２の処理で求めた乗算結果との差分を求める第３の処理とを行うためのソフトウェアを実行するソフトウェア実行部とを備えることを特徴とする多体問題計算装置。
【請求項２】
請求項１記載の多体問題計算装置において、
並行して動作する複数の力総和計算部を備え、
前記選択回路は、前記複数の力総和計算部で計算された力の総和の各座標成分を順次選択することを特徴とする多体問題計算装置。
【請求項３】
請求項２記載の多体問題計算装置において、
前記複数の力総和計算部は、計算対象粒子としてそれぞれ異なる粒子を選択し、処理単位とするイメージセルとして同一のイメージセルを選択することを特徴とする多体問題計算装置。
【請求項４】
請求項１または３記載の多体問題計算装置において、
前記力総和計算部が処理単位にしているイメージセルを示す格子ベクトルが格納される格子ベクトル保持部を備え、
前記加減算制御回路は、前記格子ベクトル保持部に保持されている格子ベクトルに基づいて、前記力総和計算部において処理単位にしているイメージセルと前記基本セルとの位置関係を判定することを特徴とする多体問題計算装置。
【請求項５】
請求項４記載の多体問題計算装置において、
前記加減算制御回路は、前記格子ベクトル保持部に保持されている格子ベクトルの各座標成分の内、前記選択回路で選択されている力の総和の座標成分と同じ座標成分が正の場合は、前記選択回路で選択された力の総和の座標成分をそのまま出力し、負の場合は符号を反転して出力し、０の場合は出力を０とすることを特徴とする多体問題計算装置。
【請求項６】
請求項１乃至５の何れか１項に記載の多体問題計算装置において、
前記第２の処理において用いる定数は、前記基本セルの辺の長さの１／２の値であることを特徴とする多体問題計算装置。
【請求項７】
多体問題計算装置の構成要素である二体間相互作用計算回路であって、
周期的境界条件を適用する系における基本セル内の粒子の１つを計算対象粒子として選択し、前記周期的境界条件のイメージセルを順次処理単位にして、そのイメージセル内の粒子が前記計算対象粒子に対して及ぼす力の総和を求めるという処理を、前記基本セル内の全ての粒子を計算対象粒子として選択するまで繰り返し行うと共に、前記基本セル内の粒子毎にその粒子に働く力の総和を求めて保持する力総和計算部と、
該力総和計算部においてイメージセル単位の力の総和が求められる毎に、その力の総和の各座標成分を順次選択する選択回路と、
前記力総和計算部において処理単位にされているイメージセルと前記基本セルとの位置関係に基づいて、前記選択回路で選択された力の総和の座標成分をそのまま出力するか、符号を反転して出力するか、或いは出力を０とする加減算制御回路と、
ビリアルの成分毎の加算結果が保持される保持回路と、
前記加減算制御回路の出力と、前記保持回路に保持されている加算結果の内の、前記加減算制御回路の出力に対応する成分の加算結果とを加算し、加算結果を前記保持回路に格納する加算器とを備えることを特徴とする二体間相互作用計算回路。
【請求項８】
請求項７記載の二体間相互作用計算回路において、
並行して動作する複数の力総和計算部を備え、
前記選択回路は、前記複数の力総和計算部で計算された力の総和の各座標成分を順次選択することを特徴とする二体間相互作用計算回路。
【請求項９】
力総和計算部が、周期的境界条件を適用する系における基本セル内の粒子の１つを計算対象粒子として選択し、前記周期的境界条件のイメージセルを順次処理単位にして、そのイメージセル内の粒子が前記計算対象粒子に対して及ぼす力の総和を求めるという処理を、前記基本セル内の全ての粒子を計算対象粒子として選択するまで繰り返し行うと共に、前記基本セル内の粒子毎にその粒子に働く力の総和を求めて保持し、
選択回路が、該力総和計算部においてイメージセル単位の力の総和が求められる毎に、その力の総和の各座標成分を順次選択し、
加減算制御回路が、前記力総和計算部において処理単位にされているイメージセルと前記基本セルとの位置関係に基づいて、前記選択回路で選択された力の総和の座標成分をそのまま出力するか、符号を反転して出力するか、或いは出力を０とし、
加算器が、前記加減算制御回路の出力と、保持回路に保持されているビリアルの各成分についての加算結果の内の、前記加減算制御回路の出力に対応する成分の加算結果とを加算し、加算結果を前記保持回路に格納し、
ソフトウェア実行部が、ビリアルの成分毎に、前記基本セル内の粒子の粒子座標と前記力総和計算部に保持されている前記基本セル内の粒子に働く力の総和とを用いて、その成分についての力の総和と粒子座標との積和を求める第１の処理と、前記保持回路に保持されているその成分についての加算結果と定数とを乗算する第２の処理と、前記第１の処理で求めた積和結果と前記第２の処理で求めた乗算結果との差分を求める第３の処理とを行うためのソフトウェアを実行することを特徴とする多体問題計算方法。
【請求項１０】
請求項９記載の多体問題計算方法において、
複数の力総和計算部が並行して動作し、
前記選択回路が、前記複数の力総和計算部で計算された力の総和の各座標成分を順次選択することを特徴とする多体問題計算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、多体問題計算技術に関する。
【背景技術】
【０００２】
分子動力学の分野では、液体、固体、高分子などの挙動を、それらを構成する原子あるいは分子の動きの結果と考え、それら粒子の動きをシミュレーションして研究が行なわれる。互いに相互作用する多数の粒子からなる系を扱う問題を多体問題といい、分子動力学では、原子や分子を粒子とした多体問題の計算を行なっている。
【０００３】
多体問題で扱う系の粒子間には相互作用が働くので、多体問題では粒子間に働く力が重要な物理量となり、必ず計算される。ある粒子に働く力は、他の全ての粒子から及ぼされる力の総和である。
【０００４】
また、物質は温度変化で圧力および体積が変化し、体積変化は粒子の位置の変化を意味するので、圧力もまた重要な物理量である。更に、その圧力を求めるために用いられるビリアルも重要な物理量として頻繁に計算される。ビリアルとは、粒子毎に、その粒子に働く力に粒子座標を乗算し、積和したものである。系内のビリアルから圧力を求めることができる。
【０００５】
ビリアルを求める装置の一例が特許文献１に開示されている。なお、特許文献１では、２つの粒子間の座標差分とその２つの粒子間に働く力とを乗算したものをビリアルとしている。
【０００６】
図８を参照すると、特許文献１に開示されている計算装置は、座標差分計算部８１と、距離二乗計算部８２と、関数計算部８３と、遅延部８４と、力計算部８５と、ビリアル計算部８６と、力総和計算部８７と、ビリアル総和計算部８８とを有している。
【０００７】
ここでは、系内に多数ある粒子の内、計算対象の粒子を粒子ｉとし、粒子ｉに影響を与える他の粒子を粒子ｊとする。また、粒子ｉの座標を（ｘi，ｙi，ｚi）とし、粒子ｊの座標を（ｘj，ｙj，ｚj）とする。粒子ｊが粒子ｉに及ぼす力の大きさをＦとする。
【０００８】
座標差分計算部８１は、２つの粒子ｉ，ｊの座標を入力とし、式（１）〜（３）に従って座標差分Δｘj，Δｙj，Δｚjを計算する。
【０００９】
Δｘj＝ｘj−ｘi  …  （１）
Δｙj＝ｙj−ｙi  …  （２）
Δｚj＝ｚj−ｚi  …  （３）
【００１０】
距離二乗計算部８２は、座標差分計算部８１で得られた座標差分から、式（４）に従って粒子間距離ｒの二乗を計算する。
【００１１】
ｒ2＝ (Δｘj)2 ＋ (Δｙj)2 ＋ (Δｚj)2   …  （４）
【００１２】
関数計算部８３は、距離二乗計算部８２で得られた粒子間距離ｒの二乗を入力とし、粒子ｊから粒子ｉに及ぼされる力Ｆを粒子間距離ｒで割った値Ｆ／ｒを計算する。粒子間に働く力は、力の種類や粒子の種類によって異なる関数形で与えられるが、いずれの場合でも粒子間距離ｒが決まれば一意に値が求まる。そのため、粒子間距離ｒの二乗から力Ｆが求まり、さらに力を距離で割った値Ｆ／ｒが求まる。
【００１３】
遅延部８４は、座標差分計算部８１で得られた座標差分Δｘj，Δｙj，Δｚjを遅延させる。この遅延は、関数計算部８３での計算にかかる時間を調整するためのものである。
【００１４】
力計算部８５は、遅延部８４からの座標差分Δｘj，Δｙj，Δｚjと、関数計算部８３からの値Ｆ／ｒとを入力とし、式（５）〜（７）に従って各座標軸方向の力の成分Ｆx，Ｆy，Ｆzを計算する。
【００１５】
Ｆx＝Ｆ×（Δｘj／ｒ）＝（Ｆ／ｒ）×Δｘj  …  （５）
Ｆy＝Ｆ×（Δｙj／ｒ）＝（Ｆ／ｒ）×Δｙj  …  （６）
Ｆz＝Ｆ×（Δｚj／ｒ）＝（Ｆ／ｒ）×Δｚj  …  （７）
【００１６】
ビリアル計算部８６は、遅延部８４からの座標差分Δｘj，Δｙj，Δｚjと、力計算部８５からの力の成分Ｆx，Ｆy，Ｆzとを入力とし、式（８）〜（１０）に従って各座標軸方向のビリアルの成分Ｖx，Ｖy，Ｖzを計算する。
【００１７】
Ｖx＝Ｆx×Δｘj  …  （８）
Ｖy＝Ｆy×Δｙj  …  （９）
Ｖz＝Ｆz×Δｚj  …  （１０）
【００１８】
力総和計算部８７は、力計算部８５からの、粒子ｉ以外の全ての粒子を粒子ｊとした、粒子ｊによる力の成分Ｆx，Ｆy，Ｆzを入力とし、全ての粒子ｊから粒子ｉに及ぼされる力の総和を求める。ビリアル総和計算部８８は、系全体のビリアルを計算する。
【００１９】
【特許文献１】特開平８−２８７０４２号公報
【発明の開示】
【発明が解決しようとする課題】
【００２０】
特許文献１に記載されている技術によれば、力およびその総和と、ビリアルおよびその総和とを求めることができる。しかし、特許文献１に記載されている技術では、力計算部において（Ｆ／ｒ）と座標差分とを乗算し、ビリアル計算部において力と座標差分とを乗算しており、力計算部だけでなくビリアル計算部にも乗算回路が必要になるため、ハードウェア量が多くなるという問題がある。このような問題点を解決するため、ビリアル計算部をソフトウェアによって実現するということが考えられる。しかし、系全体のビリアルを計算する場合、系内の粒子数Ｎに対する計算量のオーダーはＯ（Ｎ2）となり、ソフトウェアによって計算したのでは、計算時間が膨大なものになってしまう。
【００２１】
そこで、本発明の目的は、ハードウェアを用いて系全体のビリアルを計算するようにした場合には、ハードウェア量が多くなり、ソフトウェアを用いて系全体のビリアルを計算するようにした場合には、計算時間が膨大になるという課題を解決するための技術を提供することにある。
【課題を解決するための手段】
【００２２】
本発明にかかる第１の多体問題計算装置は、
周期的境界条件を適用する系における基本セル内の粒子の１つを計算対象粒子として選択し、前記周期的境界条件のイメージセルを順次処理単位にして、そのイメージセル内の粒子が前記計算対象粒子に対して及ぼす力の総和を求めるという処理を、前記基本セル内の全ての粒子を計算対象粒子として選択するまで繰り返し行うと共に、前記基本セル内の粒子毎にその粒子に働く力の総和を求めて保持する力総和計算部と、
該力総和計算部においてイメージセル単位の力の総和が求められる毎に、その力の総和の各座標成分を順次選択する選択回路と、
前記力総和計算部において処理単位にされているイメージセルと前記基本セルとの位置関係に基づいて、前記選択回路で選択された力の総和の座標成分をそのまま出力するか、符号を反転して出力するか、或いは出力を０とする加減算制御回路と、
ビリアルの成分毎の加算結果が保持される保持回路と、
前記加減算制御回路の出力と、前記保持回路に保持されている加算結果の内の、前記加減算制御回路の出力に対応する成分の加算結果とを加算し、加算結果を前記保持回路に格納する加算器と、
ビリアルの成分毎に、前記基本セル内の粒子の粒子座標と前記力総和計算部に保持されている前記基本セル内の粒子に働く力の総和とを用いて、その成分についての力の総和と粒子座標との積和を求める第１の処理と、前記保持回路に保持されているその成分についての加算結果と定数とを乗算する第２の処理と、前記第１の処理で求めた積和結果と前記第２の処理で求めた乗算結果との差分を求める第３の処理とを行うためのソフトウェアを実行するソフトウェア実行部とを備える。
【００２３】
本発明にかかる第１の二体間相互作用計算回路は、
多体問題計算装置の構成要素である二体間相互作用計算回路であって、
周期的境界条件を適用する系における基本セル内の粒子の１つを計算対象粒子として選択し、前記周期的境界条件のイメージセルを順次処理単位にして、そのイメージセル内の粒子が前記計算対象粒子に対して及ぼす力の総和を求めるという処理を、前記基本セル内の全ての粒子を計算対象粒子として選択するまで繰り返し行うと共に、前記基本セル内の粒子毎にその粒子に働く力の総和を求めて保持する力総和計算部と、
該力総和計算部においてイメージセル単位の力の総和が求められる毎に、その力の総和の各座標成分を順次選択する選択回路と、
前記力総和計算部において処理単位にされているイメージセルと前記基本セルとの位置関係に基づいて、前記選択回路で選択された力の総和の座標成分をそのまま出力するか、符号を反転して出力するか、或いは出力を０とする加減算制御回路と、
ビリアルの成分毎の加算結果が保持される保持回路と、
前記加減算制御回路の出力と、前記保持回路に保持されている加算結果の内の、前記加減算制御回路の出力に対応する成分の加算結果とを加算し、加算結果を前記保持回路に格納する加算器とを備える。
【００２４】
本発明にかかる第１の多体問題計算方法は、
力総和計算部が、周期的境界条件を適用する系における基本セル内の粒子の１つを計算対象粒子として選択し、前記周期的境界条件のイメージセルを順次処理単位にして、そのイメージセル内の粒子が前記計算対象粒子に対して及ぼす力の総和を求めるという処理を、前記基本セル内の全ての粒子を計算対象粒子として選択するまで繰り返し行うと共に、前記基本セル内の粒子毎にその粒子に働く力の総和を求めて保持し、
選択回路が、該力総和計算部においてイメージセル単位の力の総和が求められる毎に、その力の総和の各座標成分を順次選択し、
加減算制御回路が、前記力総和計算部において処理単位にされているイメージセルと前記基本セルとの位置関係に基づいて、前記選択回路で選択された力の総和の座標成分をそのまま出力するか、符号を反転して出力するか、或いは出力を０とし、
加算器が、前記加減算制御回路の出力と、保持回路に保持されているビリアルの各成分についての加算結果の内の、前記加減算制御回路の出力に対応する成分の加算結果とを加算し、加算結果を前記保持回路に格納し、
ソフトウェア実行部が、ビリアルの成分毎に、前記基本セル内の粒子の粒子座標と前記力総和計算部に保持されている前記基本セル内の粒子に働く力の総和とを用いて、その成分についての力の総和と粒子座標との積和を求める第１の処理と、前記保持回路に保持されているその成分についての加算結果と定数とを乗算する第２の処理と、前記第１の処理で求めた積和結果と前記第２の処理で求めた乗算結果との差分を求める第３の処理とを行うためのソフトウェアを実行する。
【発明の効果】
【００２５】
本発明によれば、多体問題計算装置のハードウェア量を削減し、且つ計算時間の増加を防ぐことができる。
【発明を実施するための最良の形態】
【００２６】
次に、本発明を実施するための最良の形態について図面を参照して詳細に説明する。
【００２７】
先ず、本発明の基本的な考え方について説明する。本発明は、図１に示すような格子ベクトルｎ＝（０，０，０）で示される基本セルの回りを格子ベクトルｎ＝（ｎｘ，ｎｙ，ｎｚ）で示されるイメージセルが取り囲んでいる周期的境界条件下において、ビリアルの計算式を、式（１１）に示すように系内の粒子数Ｎに対する計算量のオーダーがＯ（Ｎ２）である右辺第２項と、計算量のオーダーがＯ（Ｎ）となる右辺第１項とに分離した形式で表現する。
【００２８】
2010072680.tif 000003
【００２９】
ここで、式（１１）について詳しく説明しておく。式（１１）は、図１をｚ軸の格子ベクトルが０のＸＹ平面で見た図２において、格子ベクトルｎ＝（０，０，０）で示される各辺の長さがＬx，Ｌy，Ｌzの基本セル内のある中心粒子iの位置を座標（ｘi，ｙi,ｚi）で表し、基本セル内の粒子ｉ以外の粒子ｊの位置を座標（ｘj，ｙj，ｚj）で表し、格子ベクトルｎ＝（ｎx，ｎy，ｎz）で示されるイメージセルｎ内の粒子ｊｎを座標（ｘjn，ｙjn，ｚjn）で表し、粒子ｉに対して粒子ｊから及ぼされる力ｆijのｘ座標成分をfxijとし、粒子ｊｎから及ぼされる力ｆijnのｘ座標成分をｆxijnとし、 粒子ｉに対して全ての粒子から及ぼされる力の総和のx座標成分をＦxiとした際の、圧力テンソルのxx成分のビリアルＶxxを表している。
【００３０】
そして、本発明では、式（１１）のオーダーＯ（Ｎ）の部分の内、Ｆxiはハードウェアで計算し、Ｆxiとｘiとの乗算および乗算結果の積算はソフトウェアにより行う。また、式（１１）のオーダーＯ（Ｎ2）の部分の内のΣΣｎxΣｆxijnの部分は、イメージセルの基本セルに対する位置を示す格子ベクトルｎの値に基づいて、この格子ベクトルｎで示されるイメージセル内の力の総和を、それまでの積算値に加算あるいは減算する回路構成のハードウェアを用いて計算する。この回路構成には乗算器が不要であり、ハードウェア量の削減および消費電力の低減を図ることができる。また、オーダＯ（Ｎ2）の部分の内の定数１／２と基本セルのｘ軸方向の長さＬxとの乗算、およびこの乗算の乗算結果（１／２・Ｌx）とΣΣｎxΣｆxijnとの乗算は、ソフトウェアで計算する。
【００３１】
圧力テンソルの成分には式（１１）で示したｘｘ成分を含めて、ｙｙ、ｚｚ、ｘｙ、ｙｚ、ｚｘ、ｙｘ、ｚｙ、ｘｚの９成分があり、ｙｙ、ｚｚ成分の計算式は式（１１）の座標ｘｉ、力の総和Ｆxi、基本セルの一辺の長さＬx、格子ベクトルｎx、粒子ｊｎから及ぼされる力の和Σｆxijnをy座標成分やz座標成分に変更するだけで計算式の形式は同一である。また、ｘｙ、ｙｘ、ｙｚ、ｚｙ、ｚｘ、ｘｚ成分については、ｘｙ成分の計算式を一例として挙げると、式（１２）に示す計算式となる。
【００３２】
2010072680.tif 000004
【００３３】
式（１２）は式（１１）の座標ｘiをｙ座標ｙiに変更し、基本セルのｘ軸方向の長さＬxをｙ軸方向の長さＬyに変更し、格子ベクトルｎのｘ座標成分ｎxをy座標成分ｎyに変更したものである。ハードウェアで計算する部分ΣΣｎxΣｆxijn，ΣΣｎyΣｆxijn，ΣΣｎzΣｆxijn，…に着目すると、この部分を計算する回路は、イメージセルに含まれる粒子jによる力の各座標成分の総和Σｆxijn、Σｆyijn、Σｆzijnと、格子ベクトルの各座標成分ｎx、ｎy、ｎzとを入力し、計算する成分に応じて入力を選択する構成を採用することにより、各成分毎に計算回路を設けることなく、全ての成分を１成分分のハードウェア量の計算回路を利用して計算できることが分かる。また、力はＦxi，Ｆyi，Ｆziのように中心粒子ｉ毎に計算するのに対し、ビリアルは力と異なり系全体について総和を計算するので、パイプラインを複数用意して並列度を上げたとしても最終的に全て加算されるため、ビリアルの加算回路、加算結果の保持手段は一つで良い。以上のようにビリアルの計算回路を多体問題計算装置内に一つ持つことで、加算回路、保持回路をビリアルの全ての成分について持つ必要はなく、並列度を上げるためにパイプラインを増やしてもビリアル計算回路を増やす必要もなく、ハードウェア量の削減を可能としている。
【００３４】
〔本発明の第１の実施の形態〕
図３を参照すると、本発明の第１の実施の形態にかかる多体問題計算装置１は、二体間相互作用計算回路２と、ＣＰＵなどによって実現されるソフトウェア実行部３と、制御部４とを備えている。
【００３５】
二体間相互作用計算回路２は、中心粒子座標保持部２１と、基本セル情報保持部２２と、格子ベクトル保持部２３と、座標差分計算部２４と、座標差分遅延部２５と、距離二乗計算部２６と、関数計算部２７と、力計算部２８と、力総和計算部２９と、ビリアル計算部３０と、格子ベクトル遅延部３１とを備えている。
【００３６】
中心粒子座標保持部２１には、制御部４によって計算対象にする中心粒子ｉの座標ｒi（ｘi，ｙi，ｚi）が設定される。
【００３７】
基本セル情報保持部２２には、制御部４によって基本セルの大きさ（Ｌx，Ｌy，Ｌz）が設定される。
【００３８】
格子ベクトル保持部２３には、制御部４によって計算対象にするイメージセルを示す格子ベクトルｎ（ｎx，ｎy，ｎz）が設定される。
【００３９】
座標差分計算部２４は、制御部４から入力された計算対象にする粒子ｊの座標ｒjn（ｘjn，ｙjn，ｚjn）と、中心粒子座標保持部２１に保持されている中心粒子ｉの座標ｒin（ｘin，ｙin，ｚin）との座標差分を計算する。具体的には、周期的境界条件であるため、中心粒子ｉの座標ｒi（ｘi，ｙi，ｚi）と、粒子ｊの座標ｒjn（ｘjn，ｙjn，ｚjn）を鏡像変換した座標rjn’＝（ｘjn＋ｎxＬx，ｙjn＋ｎyＬy，ｚjn＋ｎzＬz）との差分ｒijn＝（ｘi−（ｘjn＋ｎxＬx）， ｙi−（ｙjn＋ｎyＬy），ｚi−（ｚjn＋ｎzＬz））を計算する。
【００４０】
座標差分遅延部２５は、図８に示した遅延部８４と同一の構成を有し、座標差分計算部２４で計算された座標差分ｒijnを、当該座標差分ｒijnが力計算部２８において利用されるタイミングになるまで遅延（保持）する。
【００４１】
距離二乗計算部２６は、図８に示した距離二乗計算部８２と同様の構成を有し、座標差分ｒijnから距離の二乗｜ｒijn｜2を計算する。
【００４２】
関数計算部２７は、図８に示した関数計算部８３と同様の構成を有し、ｆijn／ｒijnを計算する。
【００４３】
力計算部２８は、図８に示した力計算部８５と同様の構成を有し、ｆijn／ｒijnと座標差分ｒijnとから、粒子ｊが中心粒子ｉに及ぼす力ｆijn＝（ｆxijn，ｆyijn,ｆzijn）を計算する。
【００４４】
力総和計算部２９は、力計算部２８で計算された力ｆijnから計算対象にしている中心粒子iに及ぼされる力の総和Σｆijnをイメージセル単位で計算し、ビリアル計算部３０に出力する。また、力総和計算部２９は、計算対象にした中心粒子ｉ毎に、その粒子ｉに及ぼされる力の総和Ｆi＝（Ｆxi，Ｆyi，Ｆzi）計算し、力保持回路２９１（図４参照）に保持する。
【００４５】
格子ベクトル遅延部３１は、格子ベクトル保持部２３に保持された格子ベクトルｎを、当該格子ベクトルｎがビリアル計算部３０において利用されるタイミングになるまで遅延（保持）する。
【００４６】
ビリアル計算部３０は、力総和計算部２９で計算された力の総和Σｆijnと、格子ベクトルｎとからビリアルを計算する。
【００４７】
図４を参照すると、ビリアル計算部３０は、選択回路３０１、３０２、３０６と、加減算制御回路３０３と、加算器３０４と、ビリアル保持回路３０５とから構成されている。
【００４８】
選択回路３０１は、力総和計算部２９から入力されるセル単位の力の総和の各成分Σｆxijn，Σｆyijn，Σｆzijnを順次選択し、加減算制御回路３０３に出力する。選択回路３０２は、格子ベクトル遅延部３１から入力されている格子ベクトルｎの各座標成分ｎx，ｎy，ｎzを順次選択し、加減算制御回路３０３に出力する。
【００４９】
加減算制御回路３０３は、選択回路３０２で選択された格子ベクトルｎの座標成分の値が「正」、「負」、「０」であるかによって次のような処理を行う。
・「正」の場合；選択回路３０１で選択されたセル単位の力の総和の座標成分をそのまま加算器３０４に出力する。
・「負」の場合；選択回路３０１で選択されたセル単位の力の総和の座標成分を、符号反転して加算器３０４に出力する。
・「０」の場合；０を加算器３０４に出力する。
【００５０】
ビリアル保持回路３０５には、ビリアルの各座標成分の積算値が保持される。選択回路３０６は、ビリアル保持回路３０５に保持されているビリアルの各座標成分の積算値を順次選択し、加算器３０４に出力する。加算器３０４は、加減算制御回路３０３の出力と選択回路３０６の出力とを加算し、ビリアル保持回路３０５に保存する。
【００５１】
〔第１の実施の形態の動作の説明〕
次に、図５のフローチャートを参照して本実施の形態の動作について詳細に説明する。
【００５２】
先ず、制御部４は、基本セルの大きさＬ（Ｌx，Ｌy，Ｌz）を基本セル情報保持部２２に設定する（ステップＳ５００）。次いで、制御部４は、基本セルに含まれている粒子の中から計算対象にする粒子を１つ選択し、その粒子ｉの座標ｒi（ｘi，ｙi，ｚi）を中心粒子座標保持部２１に設定する（ステップＳ５０１）。その後、制御部４は、イメージセルの内の１つを選択し、そのイメージセルを指し示す格子ベクトルｎ（ｎx，ｎy，ｎz）を格子ベクトル保持部２３に設定する（ステップＳ５０２）。
【００５３】
その後、制御部４は、ステップＳ５０２で選択したイメージセル中の粒子を１つ選択し、選択した粒子jｎの座標rjn（ｘjn，ｙjn，ｚjn）を座標差分計算部２４に入力する（ステップＳ５０３）。
【００５４】
これにより、座標差分計算部２４において粒子ｉと粒子ｊｎとの座標差分ｒijnが計算され、力計算部２８において上記座標差分ｒijnと関数計算部２７で求められたｆijn／ｒijnとに基づいて粒子ｊｎが粒子ｉに及ぼす力ｆijnが計算され、計算結果が力総和計算部２９に出力される（ステップＳ５０４）。力総和計算部２９では、力計算部２８から出力された力ｆijnを積算する（ステップＳ５０５）。
【００５５】
制御部４は、ステップＳ５０２で選択したイメージセル内の全ての粒子ｊｎを選択するまで（ステップＳ５０６がＹとなるまで）、上記イメージセル内の粒子ｊｎの座標ｒjnを順次座標差分計算部２４に出力する。そして、上記イメージセル内の全ての粒子ｊｎを選択し、上記処理が終了すると、制御部４は、ビリアル計算部３０に対してビリアルの計算を行うことを指示する（ステップＳ５０７）。このとき、力総和計算部２９からは、ステップＳ５０２で選択したイメージセル内の粒子ｊｎが粒子ｉに及ぼす力の総和Σｆijnの各座標成分が出力されている。また、格子ベクトル遅延部３１からは、ステップＳ５０２で選択したイメージセルを指し示す格子ベクトルｎ＝（ｎx，ｎy，ｎz）が出力されている。
【００５６】
ビリアル計算部３０では、上記指示を受けると、先ず、選択回路３０２が格子ベクトルｎ＝（ｎx，ｎy，ｎz）のｘ座標成分ｎxを選択し、選択回路３０１が力の総和Σｆijnのｘ，ｙ，ｚ座標成分を順次選択し、選択回路３０６がビリアル保持回路３０５に保持されているビリアルのＸＸ，ＹＸ，ＺＸ成分を順次選択する。加減算制御回路３０３は、選択回路３０２から出力されているｎxの値が「正」であれば、選択回路３０１から順次出力される力の総和Σｆijnのｘ，ｙ，ｚ座標成分をそのまま出力し、「負」であれば選択回路３０１から順次出力される力の総和Σｆijnのｘ，ｙ，ｚ座標成分を符号反転して出力し、「０」であれば０を出力する。加算器３０４は、加減算制御回路３０３から順次出力される力の総和Σｆijnのｘ，ｙ，ｚ座標成分に関する値と、選択回路３０６から順次出力されるビリアルのＸＸ，ＹＸ，ＺＸ成分とを加算し、加算結果をビリアル保持回路３０５に再び格納する。
【００５７】
その後、選択回路３０２が格子ベクトルｎのｙ座標成分ｎyを選択し、選択回路３０１が力の総和Σｆijnのｘ，ｙ，ｚ座標成分を順次選択し、選択回路３０６がビリアルのＸＹ，ＹＹ，ＺＹ成分を選択する。これにより、加算器３０４がビリアル保持回路３０５に保持されているビリアルのＸＹ，ＹＹ，ＺＹ成分を更新する。
【００５８】
そして、最後に、選択回路３０２が格子ベクトルｎのｚ座標成分ｎzを選択し、選択回路３０１が力の総和Σｆijnのｘ，ｙ，ｚ座標成分を順次選択し、選択回路３０６がビリアルのＸＺ，ＹＺ，ＺＺ成分を順次選択する。これにより、加算器３０４がビリアル保持回路３０５に保持されているビリアルのＸＺ，ＹＺ，ＺＺ成分を更新する（ステップＳ５０８）。
【００５９】
その後、制御部４は、系内の全てのイメージセルについてビリアルの計算が完了したか否か判定する（ステップＳ５０９）。全てのイメージセルについてのビリアルの計算が完了していなければ、制御部４は、ステップＳ５０２に戻って次のイメージセルを選択する。これに対して、全てのイメージセルについてのビリアルの計算が完了していれば、基本セル内の全ての粒子ｉについてのビリアルの計算が完了しているか否か判定する（ステップＳ５１０）。
【００６０】
そして、基本セル内の全ての粒子ｉについてのビリアルの計算が完了していなければ、制御部４は、ステップＳ５０１に戻って次の粒子ｉを選択する。これに対して、全ての粒子ｉについてのビリアルの計算が完了していれば、制御部４は、二体間相互作用計算回路２（ハードウェア）を利用した計算を終了し、ソフトウェア実行部３を起動する（ステップＳ５１１）。
【００６１】
これにより、ソフトウェア実行部３は、ビリアル計算部３０内のビリアル保持回路３０５に保持されているビリアルの各成分ΣΣｎxΣｆxijn ，ΣΣｎyΣｆyijn ，ΣΣｎzΣｆzijn ，ΣΣｎyΣｆxijn ，ΣΣｎzΣｆyijn ，ΣΣｎxΣｆzijn ，ΣΣｎxΣｆyijn ，ΣΣｎyΣｆzijn ，ΣΣｎzΣｆxijnを取得すると共に、力総和計算部２９の力保持回路２９１に保持されている粒子ｉ毎の力の総和Ｆiを取得する（ステップＳ５１２）。
【００６２】
その後、ソフトウェア実行部３は、ステップＳ５１２で取得したビリアルと力の総和とに基づいて系全体のビリアルを計算する（ステップＳ５１３）。ここでは、ｘｘ、ｙｙ、ｚｚ、ｘｙ、ｙｚ、ｚｘ、ｙｘ、ｚｙ、ｘｚ成分のビリアルがそれぞれ計算される。例えば、ｘｘ成分のビリアルＶxxは、次のようにして計算される。先ず、式（１）の第１項に示す演算を行う。すなわち、粒子ｉ毎にその粒子ｉに働く力の総和のｘ座標成分とその粒子ｉのｘ座標とを乗算し、それらの和（第１の演算結果）を求める。次に、式（２）の第２項に示す演算を行う。すなわち、ビリアル計算部３０内のビリアル保持回路３０５から取得したΣΣｎxΣｆxijnに、定数１／２と、基本セルの長さＬxとを乗算し、乗算結果（第２の演算結果）を求める。その後、第１の演算結果から第２の演算結果を減算することによりビリアルＶxxを求める。
【００６３】
〔第１の実施の形態の効果〕
本実施の形態によれば、多体問題計算装置のハードウェア量を削減しつつ、計算時間の増大を防ぐことができる。その理由は、系全体のビリアルを求めるための処理を、計算量のオーダーは大きいが、乗算は不要な部分と、計算量のオーダーは小さいが、乗算が必要になる部分とに分割し、計算量のオーダーが大きくなる部分はビリアル計算部３０（ハードウェア）で計算し、計算量のオーダーが小さくなる部分はソフトウェア実行部３（ソフトウェア）で計算するようにしたからである。
【００６４】
更に、本実施の形態では、力総和計算部２９で求められた力の総和の各成分を順次選択する選択回路３０１を備えており、１つの加算器３０４を用いて複数の座標成分の計算を行うことができるので、更にハードウェア量を削減することができる。
【００６５】
〔本発明の第２の実施の形態〕
次に、本発明の第２の実施の形態について詳細に説明する。本実施の形態は、二体間相互作用計算回路を複数のパイプラインを用いて構成した場合においても、ハードウェアを削減できるようにしたことを特徴とする。
【００６６】
図６を参照すると、本発明の第２の実施の形態にかかる多体問題計算装置１ａは、二体間相互作用計算回路２、ソフトウェア実行部３および制御部４の代わりにそれぞれ二体間相互作用計算回路２ａ、ソフトウェア実行部３ａおよび制御部４ａを備えている点が、図３に示した第１の実施の形態と相違している。
【００６７】
二体間相互作用計算回路２ａは、第１および第２のパイプライン２−１、２−２と、ビリアル計算部３０ａとから構成されている。
【００６８】
第１のパイプライン２−１は、ビリアル計算部を備えていない点を除いて図３に示した二体間相互作用計算回路２と同様の構成を有している。第２のパイプライン２−２は、第１のパイプライン２−１と同様の構成を有している。
【００６９】
制御部４ａは、各パイプライン２−１、２−２内の中心粒子座標保持部２１、基本セル情報保持部２２および格子ベクトル保持部２３にそれぞれ計算対象にする中心粒子ｉの座標ｒi（ｘi，ｙi，ｚi）、基本セルの大きさＬ（Ｌx，Ｌy，Ｌz）、計算対象にするイメージセルを示す格子ベクトルｎ（ｎx，ｎy，ｎz）を設定する点、および各パイプライン２−１、２−２内の座標差分計算部２４に計算対象にする粒子ｊの座標ｒjn（ｘjn，ｙjn，ｚjn）を供給する点が、第１の実施の形態の制御部４と相違している。なお、制御部４ａは、基本セル情報保持部２２および格子ベクトル保持部２３については、各パイプラインで共通の値を設定するが、中心粒子座標保持部２１についてはパイプライン毎に異なる値を設定する。また、各パイプラインの座標差分計算部２４には同じ値を供給する。
【００７０】
ソフトウェア実行部３ａは、各パイプライン２−１、２−２内の力総和計算部２９から粒子ｉ毎の力の総和Ｆiを取得する点が、第１の実施の形態のソフトウェア実行部３と相違している。
【００７１】
図７を参照すると、本実施の形態のビリアル計算部３０ａは、選択回路３０１の代わりに選択回路３０１ａを備えている点が、図４に示した第１の実施の形態のビリアル計算部３０と相違している。
【００７２】
選択回路３０１ａは、先ず、第１のパイプライン２−１内の力総和計算部２９から入力されるセル単位の力の総和の各成分Σｆxijn，Σｆyijn，Σｆzijnを順次選択し、次いで第２のパイプライン２−２内から入力される各成分Σｆxijn，Σｆyijn，Σｆzijnを順次選択する点が、図４に示した第１の実施の形態における選択回路３０１と相違している。なお、選択回路３０１ａは、選択回路３０２が或る１つの座標成分を選択する毎に上記した処理を行う。
【００７３】
〔第２の実施の形態の動作の説明〕
次に、本実施の形態の動作について、図５のフローチャートを参照して説明する。
【００７４】
先ず、制御部４ａが、基本セルの大きさＬ（Ｌx，Ｌy，Ｌz）を各パイプライン２−１、２−２内の基本セル情報保持部２２に設定する（ステップＳ５００）。次いで、制御部４ａが、基本セルに含まれている粒子の中から、各パイプライン２−１、２−２において計算対象にする粒子をそれぞれ１つずつ選択し（パイプライン２−１用に粒子ｉ1、パイプライン２−２用に粒子ｉ2を選択したとする）、選択した各パイプライン用の粒子ｉ1，ｉ2の座標ｒi1（ｘi1，ｙi1，ｚi1），ｒi2（ｘi2，ｙi2，ｚi2）をそれぞれ各パイプライン２−１、２−２内の中心粒子座標保持部２１に設定する（ステップＳ５０１）。その後、制御部４ａは、イメージセルの内の１つを選択し、そのイメージセルを指し示す格子ベクトルｎ（ｎx，ｎy，ｎz）を、各パイプライン２−１、２−２内の格子ベクトル保持部２３に設定する（ステップＳ５０２）。
【００７５】
その後、制御部４ａは、ステップＳ５０２で選択したイメージセル中の粒子を１つ選択し、選択した粒子jｎの座標rjn（ｘjn，ｙjn，ｚjn）を各パイプライン２−１、２−２内の座標差分計算部２４に入力する（ステップＳ５０３）。
【００７６】
これにより、第１のパイプライン２−１では、力計算部２８が粒子ｊｎが粒子ｉ1に及ぼす力ｆi1jnを計算し、力総和計算部２９が上記力ｆi1jnを積算する（ステップＳ５０４、Ｓ５０５）。また、第２のパイプライン２−２では、力計算部２８が粒子ｊｎが粒子ｉ2に及ぼす力ｆi2jnを計算し、力総和計算部２９が上記力ｆi2jnを積算する（ステップＳ５０４、Ｓ５０５）。
【００７７】
制御部４ａは、ステップＳ５０２で選択したイメージセル内の全ての粒子ｊｎを選択するまで（ステップＳ５０６がＹとなるまで）、上記イメージセル内の粒子ｊｎの座標ｒjnを順次各パイプライン２−１、２−２内の座標差分計算部２４に出力する。そして、上記イメージセル内の全ての粒子ｊｎを選択すると、制御部４ａは、ビリアル計算部３０ａに対してビリアルの計算を行うことを指示する（ステップＳ５０７）。このとき、第１のパイプライン２−１内の力総和計算部２９からは、ステップＳ５０２で選択したイメージセル内の粒子ｊｎが粒子ｉ1に及ぼす力の総和Σｆi1jnの各座標成分が出力され、第２のパイプライン２−２内の力総和計算部２９からは、ステップＳ５０２で選択したイメージセル内の粒子ｊｎが粒子ｉ2に及ぼす力の総和Σｆi2jnの各座標成分が出力されている。また、格子ベクトル遅延部３１からは、ステップＳ５０２で選択したイメージセルを指し示す格子ベクトルｎ＝（ｎx，ｎy，ｎz）が出力されている。
【００７８】
ビリアル計算部３０ａでは、上記指示を受けると、先ず、選択回路３０２が格子ベクトルｎ＝（ｎx，ｎy，ｎz）のｘ座標成分ｎxを選択し、選択回路３０１ａが第１のパイプライン２−１内の力総和計算部２９から出力されている力の総和Σｆi1jnのｘ，ｙ，ｚ座標成分および第２のパイプライン２−２内の力総和計算部２９から出力されている力の総和Σｆi2jnのｘ，ｙ，ｚ座標成分を順次選択し、選択回路３０６がビリアル保持回路３０５に保持されているビリアルのＸＸ，ＹＸ，ＺＸ，ＸＸ，ＹＸ，ＺＸ成分を順次選択する。加減算制御回路３０３は、選択回路３０２から出力されているｎxの値に応じて、選択回路３０１ａから出力されている力の総和Σｆi1jnのｘ，ｙ，ｚ座標成分、力の総和Σｆi2jnのｘ，ｙ，ｚ座標成分をそのまま、或いは符号反転して、或いは値を０にして出力する。加算器３０４は、加減算制御回路３０３から順次出力される力の総和Σｆi1jnのｘ，ｙ，ｚ座標成分に関する値、力の総和Σｆi2jnのｘ，ｙ，ｚ座標成分に関する値と、選択回路３０６から順次出力されるビリアルのＸＸ，ＹＸ，ＺＸ，ＸＸ，ＹＸ，ＺＸ成分とを加算し、加算結果をビリアル保持回路３０５に再び格納する。
【００７９】
その後、選択回路３０２が格子ベクトルｎのｙ座標成分ｎyを選択し、選択回路３０１が力の総和Σｆi1jnのｘ，ｙ，ｚ座標成分、力の総和Σｆi2jnのｘ，ｙ，ｚ座標成分を順次選択し、選択回路３０６がビリアルのＸＹ，ＹＹ，ＺＹ，ＸＹ，ＹＹ，ＺＹ成分を順次選択する。これにより、加算器３０４によりビリアル保持回路３０５に保持されているビリアルのＸＹ，ＹＹ，ＺＹ成分が更新される。
【００８０】
そして、最後に、選択回路３０２が格子ベクトルｎのｚ座標成分ｎzを選択し、選択回路３０１が力の総和Σｆi1jnのｘ，ｙ，ｚ座標成分、力の総和Σｆi2jnのｘ，ｙ，ｚ座標成分を順次選択し、選択回路３０６がビリアルのＸＺ，ＹＺ，ＺＺ，ＸＺ，ＹＺ，ＺＺ成分を順次選択する。これにより、加算器３０４によりビリアル保持回路３０５に保持されているビリアルのＸＺ，ＹＺ，ＺＺ成分が更新される（ステップＳ５０８）。
【００８１】
その後、制御部４ａは、系内の全てのイメージセルについてビリアルの計算が完了したか否か判定する（ステップＳ５０９）。全てのイメージセルについてのビリアルの計算が完了していなければ、制御部４は、ステップＳ５０２に戻って次のイメージセルを選択する。これに対して、全てのイメージセルについてのビリアルの計算が完了していれば、基本セル内の全ての粒子ｉについてのビリアルの計算が完了しているか否か判定する（ステップＳ５１０）。
【００８２】
そして、基本セル内の全ての粒子ｉについてのビリアルの計算が完了していなければ、制御部４は、ステップＳ５０１に戻って各パイプライン２−１、２−２毎に次の粒子ｉを選択する。これに対して、全ての粒子ｉについてのビリアルの計算が完了していれば、制御部４は、二体間相互作用計算回路２ａ（ハードウェア）を利用した計算を終了し、ソフトウェア実行部３ａを起動する（ステップＳ５１１）。
【００８３】
これにより、ソフトウェア実行部３ａは、ビリアル計算部３０ａ内のビリアル保持回路３０５に保持されているビリアルの各成分ΣΣｎxΣｆxijn ，ΣΣｎyΣｆyijn ，ΣΣｎzΣｆzijn ，ΣΣｎyΣｆxijn ，ΣΣｎzΣｆyijn ，ΣΣｎxΣｆzijn ，ΣΣｎxΣｆyijn ，ΣΣｎyΣｆzijn ，ΣΣｎzΣｆxijnを取得すると共に、各パイプライン２−１、２−２内の力総和計算部２９から粒子ｉ毎の力の総和Ｆiを取得する（ステップＳ５１２）。
【００８４】
その後、ソフトウェア実行部３ａは、ステップＳ５１２で取得したビリアルと力の総和とに基づいて系全体のビリアルを計算する（ステップＳ５１３）。
【００８５】
〔第２の実施の形態の効果〕
本実施の形態によれば、二体間相互作用計算回路２ａを複数のパイプライン２−１、２−２を用いて構成した場合においても、ハードウェアを削減できるという効果を得ることができる。その理由は、各パイプライン２−１、２−２内の力総和計算部２９で求められた力の総和の各成分を順次選択する選択回路３０１ａを備えており、１つの加算器３０４を複数のパイプライン２−１、２−２で共用できるようにしているからである。
【産業上の利用可能性】
【００８６】
本発明は、創薬研究分野において、タンパク質の構造解析や機能解析という用途や、タンパク質と化合物の結合性の評価といった用途に利用することができる。
【図面の簡単な説明】
【００８７】
【図１】本発明が周期的境界条件をどのように利用するのかを説明するための図である。
【図２】図１の周期的境界条件をＺ軸の格子ベクトルが０のＸＹ平面で見た図である。
【図３】本発明にかかる多体問題計算装置の第１の実施の形態のブロック図である。
【図４】ビリアル計算部３０の構成例を示すブロック図である。
【図５】第１及び第２の実施の形態の動作を説明するためのフローチャートである。
【図６】本発明にかかる多体問題計算装置の第２の実施の形態のブロック図である。
【図７】ビリアル計算部３０ａの構成例を示すブロック図である。
【図８】特許文献１に記載されている計算装置を説明するためのブロック図である。
【符号の説明】
【００８８】
１、１ａ…多体問題計算装置
２１…中心粒子座標保持部
２２…基本セル情報保持部
２３…格子ベクトル保持部
２４…座標差分計算部
２５…座標差分遅延部
２６…距離二乗計算部
２７…関数計算部
２８…力計算部
２９…力総和計算部
２９１…力保持回路
３０、３０ａ…ビリアル計算部
３０１、３０１ａ、３０２、３０６…選択回路
３０３…加減算制御回路
３０４…加算器
３０５…ビリアル保持回路
３１…格子ベクトル遅延部
２、２ａ…二体間相互作用計算回路
３、３ａ…ソフトウェア実行部
４、４ａ…制御部
【図１】
2010072680.tif 000005
【図２】
2010072680.tif 000006
【図３】
2010072680.tif 000007
【図４】
2010072680.tif 000008
【図５】
2010072680.tif 000009
【図６】
2010072680.tif 000010
【図７】
2010072680.tif 000011
【図８】
2010072680.tif 000012