特許 6738087 数式処理方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6738087

(24)【登録日】2020年7月21日

(45)【発行日】2020年8月12日

(54)【発明の名称】数式処理方法

(51)【国際特許分類】

   G06F 17/10 20060101AFI20200730BHJP

【ＦＩ】

   G06F17/10 Z

【請求項の数】8

【全頁数】15

(21)【出願番号】特願2016-180966(P2016-180966)

(22)【出願日】2016年9月15日

(65)【公開番号】特開2018-45533(P2018-45533A)

(43)【公開日】2018年3月22日

【審査請求日】2019年7月4日

(73)【特許権者】

【識別番号】504261077

【氏名又は名称】大学共同利用機関法人自然科学研究機構

(74)【代理人】

【識別番号】100165663

【弁理士】

【氏名又は名称】加藤 光宏

(72)【発明者】

【氏名】伊藤 篤史

【審査官】 田中 幸雄

(56)【参考文献】

【文献】 特開平６−４４２９３（ＪＰ，Ａ）

【文献】 米国特許出願公開第２００４／０２３６８０６（ＵＳ，Ａ１）

【文献】 米国特許出願公開第２０１０／０１９１７８３（ＵＳ，Ａ１）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０６Ｆ １７／１０

(57)【特許請求の範囲】

【請求項1】

コンピュータによって数式を数学的に処理する数式処理方法であって、
前記コンピュータが実行するステップとして、
前記数式の演算子ごとの数学的な処理方法を予め記憶しておくステップと、
総和Σおよび微分の演算子を含む処理対象となる数式を入力し、前記コンピュータのメモリに記憶する入力ステップと、
前記メモリから数式を読み出し、前記記憶されている処理方法に従って、演算子に従った数学的な処理を施し、その結果を再びメモリに格納することにより演算結果としての数式を得る演算ステップとを備え、
前記演算ステップは、総和Σの変数と相関を有する変数による微分では、クロネッカーのデルタを用いて総和Σを求める式を表し、クロネッカーのデルタおよび総和Σを残した状態で演算結果としての数式を得る数式処理方法。

【請求項2】

請求項１記載の数式処理方法であって、
さらに、
前記数式の演算子ごとにコンピュータが数値計算する際の計算コードを出力するためのコード化規則を予め記憶しておくステップと、
前記演算結果としての数式に基づいて、前記コード化規則に基づいて、該数式に対応する前記計算コードを生成する計算コード生成ステップを備える数式処理方法。

【請求項3】

請求項２記載の数式処理方法であって、
前記計算コード生成ステップは、
前記数式内に含まれるクロネッカーのデルタが値１となるときに、当該デルタの係数を、該数式によって求められるべき計算結果に加算する数式処理方法。

【請求項4】

請求項２または３記載の数式処理方法であって、
前記入力ステップは、前記処理対象となる数式を複数入力可能であり、
前記演算ステップは、前記処理対象となる数式のそれぞれに対応して、前記演算結果としての数式を複数求め、
前記コード化規則は、総和Σを、計算コードの繰り返し処理であるループとして出力するための規則を含み、
前記計算コード生成ステップは、複数の前記演算結果としての数式に総和Σが共通して含まれるときには、共通のループ内で演算を行わせる計算コードを生成する数式処理方法。

【請求項5】

請求項１〜４いずれか記載の数式処理方法であって、
前記数式は、いずれも演算子によって分岐する木構造で表されている数式処理方法。

【請求項6】

請求項１〜５いずれか記載の数式処理方法であって、
前記処理対象となる数式は、分子動力学における多体ポテンシャル関数の微分を含む数式である数式処理方法。

【請求項7】

数式を数学的に処理する数式処理装置であって、
前記数式の演算子ごとの数学的な処理方法を予め記憶しておく記憶部と、
総和Σおよび微分の演算子を含む処理対象となる数式を入力する入力部と、
入力された前記数式を、前記記憶されている処理方法に従って、演算子に従った数学的な処理を施し、演算結果としての数式を得る演算部とを備え、
前記演算部は、総和Σの変数と相関を有する変数による微分では、クロネッカーのデルタを用いて総和Σを求める式を表し、クロネッカーのデルタおよび総和Σを残した状態で演算結果としての数式を得る数式処理装置。

【請求項8】

コンピュータによって数式を数学的に処理するためのコンピュータプログラムであって、
総和Σおよび微分の演算子を含む処理対象となる数式を入力し、前記コンピュータのメモリに記憶する入力機能と、
前記メモリから数式を読み出し、予め記憶された前記数式の演算子ごとの数学的な処理方法に従って、演算子に従った数学的な処理を施し、その結果を再びメモリに格納することにより演算結果としての数式を得る演算機能とをコンピュータに実現させ、
前記演算機能は、総和Σの変数と相関を有する変数による微分では、クロネッカーのデルタを用いて総和Σを求める式を表し、クロネッカーのデルタおよび総和Σを残した状態で演算結果としての数式を得る機能であるコンピュータプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、コンピュータによって数式を数学的に処理し、演算結果としての数式を得る数式処理方法に関し、特に総和Σと微分演算を含む数式の処理方法に関する。

【背景技術】

【0002】

分子動力学のシミュレーションでは、それぞれの原子に作用する力等を関数式によってモデル化することが必要となる。原子に作用する力等を表す関数は、他の多数の原子による影響を考慮した多体ポテンシャル関数の微分等によって与えられる。近年では、シミュレーションの精度向上のため、多体ポテンシャル関数がより複雑になる傾向にある。また、計算は、何百万個という多数の原子を対象として行われることになる。
分子動力学のシミュレーションを行うためには、これらの関数の微分演算等を行った数式に基づいて、数値計算するための計算コードを用意する必要があるが、多体ポテンシャルが複雑になれば、その微分を行って計算コードを用意する負荷も多大となる傾向にあった。

【0003】

数式の数学的な処理をコンピュータで行う数式処理方法については、種々の先行技術が存在する。特許文献１は、多孔質媒体内の流体のフローをモデル化するための有限要素シミュレータを作成するために、数式の記号言語翻訳機を用いてソフトウェアコードを生成する技術を開示している。また、コンピュータによる数式の処理では、Mathematica（登録商標）などの処理ソフトウェアが実用化されている。

【先行技術文献】

【特許文献】

【0004】

【特許文献1】特表２００３−５２４２２７号公報

【発明の概要】

【発明が解決しようとする課題】

【0005】

比較的簡単な多体ポテンシャルを用いて粒子に働く力を計算する例について以下に示す。ポテンシャルＵが次式で表されるとする。

【数1】

ここで、
ｒ_ｉｊは、粒子間の距離；
φは、２つの粒子間の距離によって定まる相互作用関数；
Ωは、変数ｊのとり得る範囲；
を表す。つまり、上式では、ポテンシャルＵは、ｉ番目の粒子について、Ω内にある全粒子との間でφを計算することによって求められることになる。

【0006】

かかるポテンシャルＵを考えるととき、各粒子に作用する力は、位置ｒで微分することによって求められる。この計算過程を示すと、次の通りとなる。

【数2】

【0007】

上式をさらに変形すると、次式となる。

【数3】

ここで、δは、クロネッカーのデルタであり、ｉ＝ｋまたはｊ＝ｋのときは１、その他のときは０となる関数である。

【0008】

そして、クロネッカーのデルタをさらに外すと次式となる。

【数4】

【0009】

上述の一連の計算式を図１に示した。図中の計算式（１１）がポテンシャルＵを表し、計算式（１２）〜（１６）が粒子に働く力の計算過程である。

【0010】

図２は、粒子に働く力を算出するための計算コードを例示する説明図である。図示するように、計算式（１６）は、Σで分けると２つの項からなっている。そして、第２項は、変数ｋがΩ内に存在するか否かで異なる値をとる関数となっている。
図中の枠囲みの中に、かかる計算式を計算するための計算コードを表した。図示するように、変数ｊについてのループＲ１により、計算式（１６）の第１項を計算し、変数ｉについてのループＲ２により第２項を計算する構成となっている。そして、ループＲ２においては、条件分岐を示す「if (k in Ω_ｉ) then」の文章から解る通り、変数ｋがΩに属するときにのみ計算が行われるように構成されている。

【0011】

このように、比較的単純な形式のポテンシャルを用いた計算においても、その微分を得るためには、計算式（１６）のように場合分けが必要となり、数式処理も計算コードの生成も、自動化することは困難である。分子動力学では、計算対象となるのは、力だけでなく応力など、ポテンシャルに対するさらに多様な計算も行う必要があるため、上述のような比較的単純な形式のポテンシャルであっても、要求される種々の物理量に対応する計算式を得るには、多大な時間と労力を要することになる。ポテンシャル自体の形式がより複雑になれば、こうした労力等はさらに多大となる。

【0012】

かかる課題は、分子動力学に限ったものではなく、総和Σおよび微分の演算子を含む処理対象となる数式の処理に共通の課題であった。本発明では、かかる課題に鑑み、こうした数式の処理を自動化し、数値計算のための計算コードを得るための労力の軽減を図ることを目的とする。

【課題を解決するための手段】

【0013】

本発明は、
コンピュータによって数式を数学的に処理する数式処理方法であって、
前記コンピュータが実行するステップとして、
前記数式の演算子ごとの数学的な処理方法を予め記憶しておくステップと、
総和Σおよび微分の演算子を含む処理対象となる数式を入力し、前記コンピュータのメモリに記憶する入力ステップと、
前記メモリから数式を読み出し、前記記憶されている処理方法に従って、演算子に従った数学的な処理を施し、その結果を再びメモリに格納することにより演算結果としての数式を得る演算ステップとを備え、
前記演算ステップは、総和Σの変数と相関を有する変数による微分では、クロネッカーのデルタを用いて総和Σを求める式を表し、クロネッカーのデルタおよび総和Σを残した状態で演算結果としての数式を得る数式処理方法として構成することができる。

【0014】

本発明の数式処理方法によれば、上述の通り、総和Σおよび微分の演算子をコンピュータによって数学的に処理する際に、クロネッカーのデルタを用いて総和Σを求める式を表す。こうすることによって、数式に場合分けを含むことを回避でき、統一的な形式で取扱可能となるため、数式処理を自動化することが可能となるのである。

【0015】

例えば、先に説明したポテンシャルの計算では、図１中の計算式（１５）において、以下に示すようにクロネッカーのデルタ（δ）を用いた形で数式が得られる。

【数5】

・・・（１５）
本発明では、この先を計算せず、即ち図１中の計算式（１６）の形に変形することをせず、敢えて計算式（１５）の形式でとどめるのである。図１の計算式（１５）（１６）を比較すれば明らかな通り、計算式（１５）でとどめることによって、計算式（１６）における場合分けを回避できることになる。

【0016】

計算式（１５）のクロネッカーのデルタは、ｉ＝ｋまたはｊ＝ｋのときにのみ１となり、その他は０となる関数であるから、 ０となる場合の無駄な計算をしてしまっては非効率的である。よって、従来はクロネッカーのデルタを消去した式（１６）へと変形することで、無駄な計算を行わないようにすることが一般的であった。しかし、クロネッカーのデルタの消去の結果として、前述の場合分けの発生や、総和Σの階層構造の変化が生じ、数値計算を行わせるための数式または計算コード（プログラム）を得るまでに多大な労力を要していたのである。
これに対し、本発明では、敢えてクロネッカーのデルタを残すことにより、数式処理の自動化が可能となり、数値計算のための数式を軽い負荷で得ることが可能となる。
計算式（１５）は一例に過ぎず、他の数式においても、クロネッカーのデルタを残すことにより、場合分けを回避し、総和Σの階層構造を保ったまま、統一的な形式で数式を処理できる効果は、同様に得ることができる。

【0017】

本発明においては、
さらに、
前記数式の演算子ごとにコンピュータが数値計算する際の計算コードを出力するためのコード化規則を予め記憶しておくステップと、
前記演算結果としての数式に基づいて、前記コード化規則に基づいて、該数式に対応する前記計算コードを生成する計算コード生成ステップを備えるものとしてもよい。

【0018】

こうすることにより、数値計算する際の計算コードを容易に得ることができる効果がある。また、本発明では、先に説明した通り、クロネッカーのデルタを残すことによって、場合分けを回避した統一的な形式で数式を得ることができるため、計算コードの生成も、効率的にできる利点がある。
例えば、総和Σの計算は、変数について一連の計算を繰り返すループの計算コードとして表すことができ、総和Σの階層構造はそのまま計算コード中のループの階層構造に対応する。しかし、式（１６）の様な、従来の場合分けを含む数式では、ループ内に条件分岐の計算コードを含ませる必要が生じ、全体として計算コードが複雑化する傾向にある。これに対し、本発明では、場合分けを含まない数式となるため、計算コードも概ね共通化した形式で効率的に生成することが可能となる利点がある。

【0019】

計算コードの生成は、種々の方法が可能であるが、一例として、
前記計算コード生成ステップは、
前記数式内に含まれるクロネッカーのデルタが値１となるときに、当該デルタの係数を、該数式によって求められるべき計算結果に加算するものとしてもよい。

【0020】

こうすることにより、クロネッカーのデルタは、例えば、計算式（１５）においてはｉ＝ｋまたはｊ＝ｋのときにのみ１となり、その他は０となる関数であるから、数値計算するときも、クロネッカーのデルタが値１となる場合にのみ計算を行えば足りるはずである。この点は、総和Σが含まれている場合でも同様である。上記態様では、クロネッカーのデルタが値１となる場合にのみ計算を行うようにするため、クロネッカーのデルタが０となるときに無駄な計算が行われることを回避でき、計算速度を向上させることができる。

【0021】

また、本発明においては、
前記入力ステップは、前記処理対象となる数式を複数入力可能であり、
前記演算ステップは、前記処理対象となる数式のそれぞれに対応して、前記演算結果としての数式を複数求め、
前記コード化規則は、総和Σを、計算コードの繰り返し処理であるループとして出力するための規則を含み、
前記計算コード生成ステップは、複数の前記演算結果としての数式に総和Σが共通して含まれるときには、共通のループ内で演算を行わせる計算コードを生成するものとしてもよい。

【0022】

こうすることにより、複数の数式をまとめて処理でき、より一層効率化を図ることができる。本発明では、クロネッカーのデルタを含む統一的な形式で数式を処理し、数式ごとに個別の複雑な処理を回避できるからこそ、このように複数の数式をまとめて処理できるのである。
即ち、従来の方法では、クロネッカーのデルタの消去において総和Σの階層構造が変わってしまうために、対象となるポテンシャルと、力などの各種微分量の間で、総和Σの階層構造が異なってしまう。これに対し、本発明の方法では、クロネッカーのデルタの消去を行わないために、ポテンシャルとその各種微分量の間で、総和Σの階層構造が同一のまま保たれる。よって、計算コードにおいてループの融合ができ、効率的な計算が可能となる。

【0023】

本発明は、種々の形式で表した数式を対象とすることができるが、特に、
前記数式は、いずれも演算子によって分岐する木構造で表されているものとしてもよい。

【0024】

こうすることによって、再帰処理によって数式を処理することが可能となり、本発明をコンピュータによって実現するためのプログラム構造の簡略化を図ることが可能となる。

【0025】

本発明は、種々の分野の数式の処理に適用可能であるが、
特に、前記処理対象となる数式は、分子動力学における多体ポテンシャル関数の微分を含む数式とすることが有用である。

【0026】

分子動力学では、多体ポテンシャル関数を微分等して種々の物理量を演算することになり、また多体ポテンシャル関数がより複雑になる傾向にある。しかも、計算は、何百万個という多数の原子を対象として行う必要がある。かかる計算に要する負荷は非常に大きなものとなるため、本発明により、特に数式処理の負荷を軽減できる効果は非常に大きい。従って、本発明は、特に分子動力学の分野に有用性が高い。

【0027】

本発明について、上述した種々の特徴は、必ずしも全てを備えている必要はなく、適宜、その一部を省略したり組み合わせたりしてもよい。
また、本発明は、上述した数式処理方法で説明した種々の機能をコンピュータによって実現するためのコンピュータプログラムとして構成してもよいし、かかるコンピュータプログラムを記録したコンピュータ読み取り可能な記録媒体として構成することもできる。ここで記録媒体としては、ハードディスク、メモリ、ＣＤ−ＲＯＭ、フラッシュメモリなどの媒体を利用することができる。
本発明は、さらに、上述した数式処理方法を実現する数式処理装置として構成することもできる。かかる装置は、例えば、数式処理方法で説明した種々の機能を実現するためのコンピュータプログラムをコンピュータにインストールすることによってソフトウェア的に構成することもできるし、各機能を専用のハードウェア回路によって構成することもできる。

【図面の簡単な説明】

【0028】

【図1】ポテンシャルから粒子に作用する力を得る一連の計算式を示す説明図である。

【図2】粒子に働く力を算出するための計算コードを例示する説明図である。

【図3】自動コード生成処理のフローチャートである。

【図4】関数例を示す説明図である。

【図5】数式を木構造で表した例を示す説明図である。

【図6】計算コードの出力例（１）を示す説明図である。

【図7】計算コードの出力例（２）を示す説明図である。

【図8】計算コードの出力例（３）を示す説明図である。

【図9】計算コードの評価結果を示す説明図である。

【発明を実施するための形態】

【0029】

以下、コンピュータによって数式を処理し、数値計算用の計算コードを自動で得る数式処理装置としての例に基づき、本発明の実施例について説明する。
本実施例では、数式処理装置は、以下で説明する種々の機能を実現するためのコンピュータプログラムを、コンピュータにインストールすることによってソフトウェア的に構成している。機能の少なくとも一部を、専用の回路などのハードウェアによって構成する方法をとってもよい。また、スタンドアロンで稼働する一台のコンピュータで構成してもよいし、ネットワークで接続された複数のコンピュータによって構成してもよい。

【0030】

Ａ．関数の読み込み：
図３は、自動コード生成処理のフローチャートである。実施例の数式処理装置が行う全体処理のフローチャートに相当する。
処理を開始すると、コンピュータは、まず処理対象となる関数を読み込む（ステップＳ１０）。図の右側に関数の例を示した。分子動力学上の処理を行う場合には、ポテンシャルエネルギーのモデルとともに、それを微分して得られる粒子に働く力、力と位置ベクトルとの内積で得られる応力などが処理対象となる関数として挙げられる。
ここで、「処理」とは、数式を数学的に変形することを意味する。図中の例で言えば、粒子に働く力として、ポテンシャルエネルギーを位置で微分する数式が入力されると、この微分演算を数学的に実行し、計算結果としての数式を得ることを言う。

【0031】

図中では、通常の数学的表記で数式を表しているが、それぞれの関数は、コンピュータが入力可能な形式で表記すればよい。かかる表記としては、例えば、ＬａＴｅＸやＡＭＳ−ＬａＴｅＸなどを利用して、各数式をテキスト形式で表記する方法をとることができる。

【0032】

図４は、関数例を示す説明図である。分子動力学で利用される関数の例を示した。もっとも、本実施例における処理対象となる関数は、これらに限定されるものではない。
図の左上には、全ての関数の基礎となるポテンシャルエネルギー（関数１）を示した。ここでは、Ｕと表記してあるが、実際には、粒子の位置座標や種々のパラメータａを含む、具体的な関数として定義されることになる。
ポテンシャルエネルギー（関数１）については、そのパラメータａによる微分（関数２）の計算が求められることもある。
また、ポテンシャルエネルギー（関数１）を座標によって微分することにより、粒子に働く力（関数３）が得られる。
この粒子に働く力（関数３）については、任意のベクトルとのテンソル積（関数４）や、そのパラメータａによる微分（関数５）の計算が求められることもある。
さらに、粒子に働く力（関数３）と位置とのテンソル積によって応力（関数６）が得られ、そのパラメータａによる微分（関数７）の計算が求められることもある。
本実施例では、これらのいずれの関数も処理対象とすることができる。即ち、図３のステップＳ１０において、図４に示した各関数を入力すれば、ポテンシャルエネルギーの内容に基づいて数学的な演算が行われることになる。

【0033】

本実施例では、数式を木構造で取り扱う。木構造について、説明する。
図５は数式を木構造で表した例を示す説明図である。木構造では、図示するように数式をまとまりごとにノードとして捉える表記方法である。
図の例は、最上位の関数、即ち数式全体をφとし、その内容が、「Ａ＋Ｂ＋Ｃ」で表されることを意味している。φおよび「Ａ＋Ｂ＋Ｃ」がそれぞれノードとなる。
次に、「Ａ＋Ｂ＋Ｃ」の各項ごとにみると、図の例では、Ａの内容が、「ａ×ｂ×ｃ」であることを表している。この「ａ×ｂ×ｃ」もノードとなる。図示を省略しているが、同様に、Ｂ、Ｃも何らかの数式となる。Ｂ，Ｃにも、それぞれ対応するノードが存在することになる。
さらに、「ａ×ｂ×ｃ」については、ａが「３」という数値、ｂが「ｆ（ｒ）」という関数、ｃが「ｒ」という変数で表され、それぞれがノードを構成する。ｂに対応する関数「ｆ（ｒ）」については、その内容に応じて、さらに下位のノードが存在することになる。
図中の例では、結局、「ａ×ｂ×ｃ」は、「３×ｆ（ｒ）×ｒ」という数式を表すことになり、「Ａ＋Ｂ＋Ｃ」およびφは、「３×ｆ（ｒ）×ｒ＋Ｂ＋Ｃ」という数式を表すことになる。

【0034】

このように木構造で数式を表現すると、それぞれの定型的な演算を処理するモジュールを用意し、これを再帰的に呼び出すことで、全体の数式処理を比較的簡易なプログラム構造で実現することが可能となる。

【0035】

例えば、図の例において、φを微分することを考える。この演算は、ノード「Ａ＋Ｂ＋Ｃ」という多項式の微分演算を行うモジュールを呼び出すことになり、それぞれＡ，Ｂ，Ｃの各項を微分するという処理結果となる。ただし、Ａ，Ｂ，Ｃの微分演算自体はこの時点では完了していない。
そこで、上記結果を受け、今度は、下位のノードであるＡ、即ち「ａ×ｂ×ｃ」の微分を行うモジュールを呼び出す。その後は、結果に応じて、さらに、ａ，ｂ，ｃの各ノードの微分を行うモジュールを呼び出すことになる。このように演算が完了するまで、下位のノードに計算を進め、演算のモジュールを再帰的に呼び出すことになるのである。
こうして、下位のモジュールを呼び出す必要がないノード（例えば、図中のａのノードなど）に到達すると、そのノードについての演算結果が確定する。すると、その結果を上位のノードに引き渡すことによって、順次、演算結果が確定し、最終的には、最上位のノードであるφの微分結果が得られるのである。
木構造を採用すると、このように再帰的呼出により、数式処理を比較的容易に行うことが可能となる利点がある。

【0036】

Ｂ．微分演算の実行：
図３に戻り、自動コード生成処理について説明する。コンピュータは、関数の読み込み（ステップＳ１０）を終えると、入力された関数に基づいて微分演算を実行する（ステップＳ２０）。図の右側に、粒子に働く力の微分演算例を示した。計算過程は、図１に示した通りである。本実施例では、こうした関数の数学的な処理をコンピュータによって実現するのである。

【0037】

本実施例では、図５に示したように、数式を木構造で捉え、ノードごとに演算を実行するようにした。具体的には、ノードの種類ごとに微分演算を実行するためのモジュールを用意し、これを再帰的に呼び出すことによって演算を実行する。本実施例で用意されているモジュールを以下に例示する。

【0038】

（１）数値型のノード：
数値型のノードとは、図５中のノードａのように数値のみからなるノードである。かかるノードの微分演算に対しては、常に微分不能を意味するＮＵＬＬを結果として返還するモジュールを用意した。

【0039】

（２）スカラー変数型のノード：
スカラー変数型のノードとは、図５中のノードｃのようにスカラーの変数からなるノードである。かかるノードの微分演算に対しては、次のように変数の定義に応じて演算結果を返還するモジュールを用意した。
まず、ノードの変数名が微分変数と同じ場合は微分可能と判断する。この場合、ノードの変数名がグローバルパラメータの場合はクロネッカーデルタ型ノードを返し、添え字がある場合はクロネッカーデルタに添え字をつける。また、ユーザー定義関数の引数の場合は単に数値型ノードの1を返す。
一方、ノードの変数名が微分変数と異なる場合は、微分不可能を表すＮＵＬＬを返す。

【0040】

図１に示した数式では、式（１４）中のカッコ内のｒ_ｉ，ｒ_ｊをそれぞれｒ_ｋで微分する演算が、このスカラー変数のノードに対する演算に相当する。従って、図１中の式（１５）のように、クロネッカーのデルタが生成されているのである。

【0041】

（３）線形結合型のノード：
線形結合型のノードとは、図５中の「Ａ＋Ｂ＋Ｃ」のように関数の和の形で表されるノードである。かかる場合の微分演算に対しては、Ａ，Ｂ，Ｃという下位のノードをそれぞれ微分して線形結合した結果を返すモジュールを用意した。下位のノードの微分は、その内容に応じたモジュールが再帰的に呼び出されることになる。また、下位のノードの微分結果が全てＮＵＬＬの場合には、線形結合型のノードの微分結果もＮＵＬＬとなる。

【0042】

（４）積結合型のノード：
積結合型のノードとは、図５中の「ａ×ｂ×ｃ」のように関数の積の形で表されるノードである。かかる場合の微分演算に対しては、（ａ×ｂ×ｃ）’＝ａ’×ｂ×ｃ＋ａ×ｂ’×ｃ＋ａ×ｂ×ｃ’（ただし、’は微分を表す）というように処理するモジュールを用意した。この処理には、ａ，ｂ，ｃという下位のノードの微分が必要となるため、それぞれ下位のノードの内容に応じたモジュールが再帰的に呼び出されることになる。

【0043】

同様にして、種々のノードごとに微分演算を実行するモジュールを用意する。ノードの種類としては、例えば、ベクトル内積を表すノード、ベクトル外積を表すノード、総和（Σ）を表すノード、三角関数や指数関数などの各種数学関数を表すノードなどが挙げられる。

【0044】

本実施例では、以上の演算において、クロネッカーのデルタが生成されたとき、クロネッカーのデルタに対しては、それ以上の演算は行わない。図１の一連の計算式について説明すれば、計算式（１５）において、クロネッカーのデルタを演算して計算式（１６）の形への変形は行わず、計算式（１５）の形で演算を停止する。

【0045】

図１の計算式（１５）（１６）を比較してわかる通り、計算式（１５）には総和Σが二重で作用しているのに対し、クロネッカーのデルタを外す場合には、計算式（１６）に示されるように総和Σを一つ外すことができる。かかる意味では、クロネッカーのデルタを外せば、式を簡略化することができるのであるが、本実施例では、敢えてこれを外さずに、二重のΣが付された計算式（１５）の状態で演算を停止するのである。

【0046】

Ｃ．計算コードの出力：
図３に戻り自動コード生成処理について説明する。
コンピュータは、ステップＳ２０の演算処理を実行すると、次に数値計算用プログラムコード（これを計算コードと呼ぶこともある）の出力を行う（ステップＳ３０）。ステップＳ２０での演算は、数学的な演算処理であるのに対し、ステップＳ３０では、演算処理の結果を用いて、数値計算するための計算コードを生成するのである。
図の右側に計算コードの出力例を示した。このように、ステップＳ３０は、プログラムコードがアウトプットとなる。

【0047】

図６は、計算コードの出力例（１）を示す説明図である。図１中の計算式（１５）に対する計算コードの例を示した。計算式（１５）はｋ番目の粒子に働く力を求める式である。計算コードでは、全ての取り得るｋ＝1〜N、（Nは粒子数）、について、粒子に働く力を計算する。計算の実行時に、ｋが引数として与えられるものとして計算コードを生成してもよい。計算コードでは、初期化処理として、ポテンシャル、力Ｆのバッファをクリアした後、計算式（１５）の計算に入るようになっている。

【0048】

図示するように、計算式（１５）には、ｉ，ｊそれぞれについて総和（Σ）が二重に作用している。そこで、計算コードでは、これをｊによる繰り返しループ、ｉによる繰り返しループの二重ループに対応づけるものとした。
また、全てのｋについて計算を行うため、ｉ，ｊの外側に、ｋによる繰り返しループが設けられている。

【0049】

ループ内の処理Ｃ１の部分では、図中の変数ｔｅｍｐに対応する式により、それぞれｉ，ｊの値に応じて、クロネッカーのデルタの係数に相当する部分を計算する。
そして、クロネッカーのデルタの処理を行う。δ_ikが意味を持つのは、ｉ＝ｋのときのみであるから、「if(i==k)」つまり「ｉ＝ｋが真」という条件が成立するときに、力Ｆ［ｋ］のバッファにｔｅｍｐの計算結果を加えればよい。また、同様に−δ_jkが意味を持つのは、ｊ＝ｋのときのみであるから、「if(j==k)」つまり「j＝ｋが真」という条件が成立するときに、力Ｆ［ｋ］のバッファに−ｔｅｍｐの計算結果を加えればよい。

【0050】

図６に示すように計算コードを生成してもよいが、上述の通り条件分岐（ｉｆ文）が含まれるため、計算コードを更に簡略化することを考える。
図７は、計算コードの出力例（２）を示す説明図である。この計算コードでは、主たる計算Ｃ２において、図６の計算コードに対して、変数ｋについてのループが、ｉ，ｊのループの下位に入れ換えられている。変数ｋは、ｉ，ｊと独立の変数であるから、このようにループの順序を入れ換えても結果は変わらない。
この計算Ｃ２において、ｉ＝ｋなる条件分岐（ｉｆ（ｉ＝＝ｋ））では、結局、Ｆ［ｋ］＝Ｆ［ｋ］＋ｔｅｍｐなる計算が行われるのは、ｉ＝ｋのときのみである。従って、条件分岐を外し、Ｆ［ｉ］＝Ｆ［ｉ］＋ｔｅｍｐと記載しても計算結果は変わらない。同様に、ｊ＝ｋなる条件分岐（ｉｆ（ｊ＝＝ｋ））についても、ｋをｊに置換して、Ｆ［ｊ］＝Ｆ［ｊ］＋ｔｅｍｐと記載しても計算結果は変わらない。ただし、ｋについてのループを残しておくと、ｉ＝ｋなる条件分岐（ｉｆ（ｉ＝＝ｋ））およびｊ＝ｋなる条件分岐（ｉｆ（ｊ＝＝ｋ））の判定処理が、ｋのループに従って繰り返し行われることになるから、上述の置換は、ｋのループを外した上で行う必要がある。
こうして変形された計算コードの出力例が、図８に示す計算コードの出力例（３）である。ここでは、ｋについてのループは外れている。また、ｉ＝ｋなる条件分岐（ｉｆ（ｉ＝＝ｋ））およびｊ＝ｋなる条件分岐（ｉｆ（ｊ＝＝ｋ））も外れている。このように、全てのｋについて計算を行うときは、ｋについてのループを設けるまでなく、結果を得ることができるのである。このとき、計算式（１５）におけるクロネッカーのデルタは、それぞれＦ［ｉ］、Ｆ［ｊ］に係数の計算結果ｔｅｍｐまたは−ｔｅｍｐを加える処理として実現されることになる。

【0051】

このように総和Σを繰り返しループに対応させ、クロネッカーのデルタを計算結果のバッファへの足し込みに対応させることにより、計算コードを生成することができる。それぞれｉ，ｊが必要な粒子の範囲をループによって移り変わる為、全ての取り得るｋ＝1〜Nについての計算をしたことと等価であることが数学的に保障される。
図６〜図８では、計算式（１５）を例にとって式に忠実な形の計算コード例（図６）から、簡略化された計算コード例（図８）を得る過程を示したが、クロネッカーのデルタの処理は、計算式（１５）に限らず同様に成立するものであるから、本実施例では、図６の形、即ちｋについての繰り返しループを設ける計算コードではなく、図８の形、即ちｋについての繰り返しループを外した計算コードを出力するようにしている。
本実施例では、関数を木構造で表しており、その演算結果も木構造で捉えることができる。従って、繰り返しループ内での計算コードの内容は、木構造に従って、順次、下位のノードを呼び出すことにより生成することができる。

【0052】

図６には、計算式（１５）に対応する計算コードを例示した。本実施例では、先に図４に示したように多様な関数を処理することができる。また、分子動力学では、これらの関数の演算結果は、総和Σが三重以上に作用する多重構造で表されることが多い。特に、本実施例では、クロネッカーのデルタを外さない状態で演算をとめるため、このような統一的な形式で演算結果を得ることができる。
従って、本実施例において、複数の関数を処理する場合には、図６中の処理Ｃ１の部分に、各関数の計算内容を含めるようにしてもよい。例えば、入力された複数の関数の演算結果を比較し、総和Σが同じ形式で含まれているもののみを抽出して、処理Ｃ１の部分に含めるようにしてもよい。また別の方法として、一旦、関数ごとに計算コードを生成した後、同様の繰り返しループとなる部分をまとめるようにしてもよい。
このように一つの繰り返しループ内で複数の関数を計算可能とすることにより、ｉ、ｊに対する繰り返しループ自体を、関数の数だけ重複して行うことを回避でき、より効率的に数値計算可能な計算コードを得ることが可能となる。

【0053】

図９は、計算コードの評価結果を示す説明図である。ＥＡＭ（Embedded Atom Method）と呼ばれるポテンシャルモデルとして、タングステンーヘリウム系の相互作用を記述する多体ポテンシャルを用いて評価した結果である。三角、丸、四角は、それぞれ粒子数が異なる結果を示しており、この順に、粒子数が増大している。また、横軸は、この計算を行うコンピュータのＣＰＵのコア数であり、縦軸は単位時間内に計算されるステップ数である。つまり、縦軸は上ほど、処理が速いことを表している。
図示する通り、いずれの粒子数においても、本実施例による自動生成による処理速度は、オペレータがマニュアルで計算コードを生成した場合の処理速度に匹敵する結果となっていることが解る。従って、本実施例による自動生成は、十分に実用的な計算コードを生成していることが示されている。

【0054】

以上で説明した本実施例によれば、総和Σおよび微分の演算子を含む関数の数学的な演算を効率的に行うことができる。また、クロネッカーのデルタが表れた時点で、クロネッカーのデルタについては、それ以上の演算を行わないことにより、対象となるポテンシャル関数とその各種微分量で演算結果のループ構造を統一的な形式に整えることができ、これに基づく計算コードを自動生成することが可能となる。
この結果、分子動力学のように多数の粒子を対象とする複雑な計算式をコンピュータで数値可能な状況を整えるための労力を大幅に軽減することができる利点がある。

【0055】

以上の説明では、具体的な関数を例示したが、本実施例は、多種多様な関数について適用可能である。また、処理対象となる関数は、分子動力学に限られない。

【産業上の利用可能性】

【0056】

本発明は、コンピュータによって、総和Σと微分演算を含む数式を処理し、数値計算用の計算コードを得るために利用することができる。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】