(57)【要約】【課題】分子軌道法におけるクーロン力の様に空間内の位置によって値が著しく変化する多変数関数を数値積分する際に、必要な精度を得る事が困難だった。【解決手段】空間Aにおいてほぼ一様に分布する標本点ｐから被積分関数が定義された空間Bの標本点ｑに対応付けを行う際し、点ｐに対する空間要素ｖと点ｑに対する空間要素ｗの比を改良して、分子軌道法における原子核の位置の様に空間Bにおける被積分関数の値ｆ（ｑ）が特に大きくなる点ｑ０の近傍の標本点ｑでも数値積分の項ｆ（ｑ）ｗが特に大きくならない様にしたため、特定の標本点のｆ（ｑ）ｗの誤差が積分値全体に対して支配的になるのを防ぐことが出来た。その結果、GeSe分子の様な原子番号が比較的大きな原子を含む分子においても、各原子に働く力を図２に示す様に正しく計算できるようになった。【選択図】図２
2014182418.tif 000002
【特許請求の範囲】
【請求項１】
空間Bで定義された被積分関数の特異点ｑｊ，０毎に対応する空間Aｊ（ｊ＝１，２，・・・，Ｍ）を設定し、
各空間Aｊにおいて、
所定の方式でNｊ個の標本点ｐｊ，ｉ（ｉ＝１，２，・・・，Ｎｊ）を配列させる工程と、
標本点ｐｊ，ｉに対して重みｖｊ，ｉを与える工程と、
空間Aｊから空間Bへの標本点の変換式を設定する工程と、
標本点ｐｊ，ｉを、該変換式により変換しさらに特異点ｑｊ，０の位置だけ移動して空間Bの標本点ｑｊ，ｉに対応付ける工程と、
標本点ｑｊ，ｉにおける空間Aと空間Bの標本点の重みの比ｄｗ/ｄｖと被積分関数の空間Aｊへの分配比ｇｊを算出し、
標本点ｑｊ，ｉの重みＷｊ，ｉとして
【数２４】
2014182418.tif 000029
を算出する工程を有する多変数関数の数値積分用の標本点とその重みの生成方法において、
該重みの比ｄｗ/ｄｖは、特異点ｑｊ，０と標本点ｑｊ，ｉとの距離ｒに対してｒ２ｈ（ｒ）の形で依存し、
関数ｈ（ｒ）はｒ→０に対して減少しつつ所定の値に収束する関数であり、
該所定の値はｈ（ｒ）が含むパラメーターの調整によって０または任意の正の値に出来るものである事を特徴とする、
多変数関数の数値積分用の標本点とその重みの生成方法。
【請求項２】
請求項１に記載の多変数関数の数値積分用の標本点とその重みの生成方法において、
前記関数ｈ（ｒ）は、０又は正の値を取るパラメーターR０を含む
【数２５】
2014182418.tif 000030
である事を特徴とする多変数関数の数値積分用の標本点とその重みの生成方法。
【請求項３】
空間Bで定義された被積分関数のM個の特異点ｑｊ，０の座標を受け取り、
各空間Aｊ（ｊ＝１，２，・・・，M）において
空間Aｊから空間Bへの標本点の変換式を設定し、
標本点の個数Ｎｊを設定し
所定の方式に従って標本点ｐｊ，ｉ（ｉ＝１，２，・・・，Ｎｊ）を生成し、
標本点ｐｊ，ｉに付随する重みｖｊ，ｉを設定し、
前記変換式に従って標本点ｐｊ，ｉを変換し特異点ｑｊ，０の位置だけ移動し空間Bの標本点ｑｊ，ｉに対応付し、
標本点ｑｊ，ｉにおける空間Aと空間Bの標本点の重みの比ｄｗ/ｄｖと被積分関数の空間Aｊへの分配比ｇｊを算出し、
標本点ｑｊ，ｉに付随する重みｗｊ，ｉとして数式２６を算出し、
【数２６】
2014182418.tif 000031
空間Bの標本点ｑｊ，ｉとその重みｗｊ，ｉを出力する多変数関数の数値積分用の標本点とその重みの生成プログラムにおいて、
該重みの比ｄｗ/ｄｖは、空間Bの特異点ｑｊ，０と標本点ｑｊ，ｉとの距離ｒに対してｒ２ｈ（ｒ）の形で依存し、
関数ｈ（ｒ）はｒ→０に対して減少しつつ所定の値に収束する関数であり、
該所定の値はｈ（ｒ）が含むパラメーターの調整によって０または任意の正の値に出来るものである事を特徴とする、
多変数関数の数値積分用の標本点とその重みの生成プログラム。
【請求項４】
請求項３に記載の多変数関数の数値積分用の標本点とその重みの生成プログラムにおいて、
前記関数ｈ（ｒ）は、０又は正の値を取るパラメーターR０を含む
【数２７】
2014182418.tif 000032
である事を特徴とする多変数関数の数値積分用の標本点とその重みの生成プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、数値シミュレーション等において必要とされる多変数関数の数値積分を高精度で行うために好適な標本点と重みの生成方法及びその生成のためのプログラムに関するものである。
【背景技術】
【０００２】
数値シミュレーション等において、必要とされる関数の積分が解析的には不可能または困難な場合がある。その場合でも被積分関数が１変数の関数であれば、Ｓｙｍｐｓｏｎの公式を始めとする各種公式を用いて任意の精度で数値積分を行う事が出来る。しかし被積分関数が多変数関数の場合は、単純に重積分を行ったのでは計算量が著しく増大するので、モンテカルロ法が用いられる事が多い。モンテカルロ法とは、非特許文献１に記載されている様に空間の点ｐにおける被積分関数をf（ｐ）とする時、数式１に示す積分を数式２に示す形で近似する方法である。
【０００３】
【数１】
2014182418.tif 000003
【０００４】
【数２】
2014182418.tif 000004
【０００５】
ここでpiは積分を行う空間に分布するN個の標本点の一つで、例えばｆ（p）が３変数の関数ならｐｉ＝（ｓｉ，ｔｉ，ｕｉ）の様に３個の座標値を持つ。またvｉは標本点piに割り当てられる重み（空間要素の体積）である。標本点としては準乱数の組が広く用いられている。準乱数とは、適当な無理数を整数倍した値の小数点以下の部分であり、整数として１〜Nを順次与える事によってN個の準乱数が得られる。また３変数の関数の場合は、例えば√２、√３、√５の様な相異なる３個の無理数に対し整数として１〜Nを順次与える事によって標本点がN個得られる。この様にして得られた標本点は空間の中で一様に近い分布をする。通常は完全に一様な分布と見なしてｖｉ＝Ｖ／Ｎを与える。ここでVは積分を行う空間の体積である。数式２の誤差は概ねNに反比例するので、高い精度が必要なら通常多数の標本点を使う必要がある。
【０００６】
数値積分が必要となる関数の中には、積分を行う空間の中でその値（絶対値）を大きく変えるものがある。例えば数値積分が必要となる典型的な数値シミュレーションである分子軌道法の場合、被積分関数として数式３の形の関数が現れる事が多い。
【０００７】
【数３】
2014182418.tif 000005
【０００８】
ここでＺは分子を構成する原子の原子番号、ｒは原子核と標本点との距離、ｎは通常１か１以上である。例えば分子の全エネルギーを求めるには原子核と電子のクーロン相互作用のエネルギーを計算する必要があるが、その場合はｎ＝１である。１/ｒｎの寄与に加えて、波動関数φ１及びφ２は原子核に近づくにつれ値が大きくなるので数式３は原子核近傍で著しく大きな値を取る事が多い。原子核の位置の様に被積分関数の値が著しく大きくなる点を以下では被積分関数の特異点、又は単に特異点と呼ぶ。一方原子核から遠くなると一般に波動関数は０に向かって指数関数的に収束する。そのため数式２の右辺の各項の大きさは著しく不揃いとなり、大きな値をとる項の誤差が支配的となってNを大きくしてもその割には積分の精度が向上しない。この効果は特にZの大きな原子を含む分子の場合に顕著となる。
【０００９】
この問題を改善するため変数の変換が利用される事が多い。始めに原子が1個だけある場合について説明する。始めに仮想的な空間Aを考え、この空間における標本点piに対し、適当な変換式によって被積分関数が定義された空間Bにおける標本点ｑｉを対応させる。所望の積分は空間Bにおいて数式４の形の数値積分で近似される。
【００１０】
【数４】
2014182418.tif 000006
【００１１】
ｗiは標本点ｑｉの重みである。空間Aと空間Bの標本点の重みの比ｄｗ/ｄｖを使うとｗiは数式５の様に表される。
【００１２】
【数５】
2014182418.tif 000007
【００１３】
ここで変換式を調整して、空間Bの標本点を、関数ｆ（ｑ）の値が大きい場所では密に分布させ、ｆ（ｑ）の値が小さい場所では疎に分布させることにする。するとｆ（ｐ）の値が大きい場所ではｄｗ/ｄｖは小さくなり、ｆ（ｐ）の値が小さい場所ではｄｗ/ｄｖが大きくなるため、数式４の右辺の各項の値のバラツキが小さくなる。このため数式４を直接数値積分するより、Nが同じでも精度が高くなる。見方を変えると、空間Bにおける被積分関数ｆ（ｐ）を、空間Aにおいて値の変化が少ない関数ｆ（ｐ）ｄｗ/ｄｖに変換して数値積分すると考える事も出来る。
【００１４】
３変数の場合について具体的な例を考えてみる。非特許文献２にならって空間Aとしてｓ＝[０，１]、ｔ＝[０，１]、ｕ＝[０，１]の立方体を想定し、この中に準乱数を使って標本点を分布させる。空間Bの座標系として極座標系を採用し、数式６の変換式に従って標本点pi（ｓｉ，ｔｉ，ｕｉ）を標本点ｑi（Φｉ，θｉ，ｒｉ）に対応付ける。
【００１５】
【数６】
2014182418.tif 000008
【００１６】
ここでR０は標本点の分布の広がりを調整するパラメーターである。Γはｕ→１とする時ｒ→∞となるように定められる係数であり数式７で与えられる。この時ｄｗ/ｄｖは数式８の様に与えられる。
【００１７】
【数７】
2014182418.tif 000009
【００１８】
【数８】
2014182418.tif 000010
【００１９】
従って関数の特異点（原子核の位置）の近傍で被積分関数ｆ（ｑ）が数式３に従って大きくなっても、ｎ＝１であれば数式４の右辺各項は特に大きくならないので数値積分の精度は低下し難い。
【００２０】
以下ではM個の原子を含む一般の分子について説明する。この場合、空間BにM個の特異点ｑｊ，０（原子核の位置，ｊ＝１，２，・・・，Ｍ）を含む事になる。空間Aの標本点からM個の特異点を含む空間Bの標本点に、単一の変換式で対応づけられると理想的であるが、今のところ適当な変換式が見つかっていない。そこで各特異点に対応する空間A１、空間A２、・・・、空間AMを考え、各空間において所定の方式に従って標本点ｐｊ，ｉ（ｊ＝１，２，・・・，Ｍ，ｉ＝１，２，・・・，Ｎｊ）を配置することにする。Ｎｊは特異点ｑｊ，０の性質によって最適な値が異なるので一般にはｊに依存する。次に数式６に従って標本点ｐｊ，ｉから空間Bに変換する。各空間によってR０や係数Γも異なっていて良い。さらに特異点ｑｊ，０の座標分だけ平行移動し標本点ｑｊ，ｉとする。
ここで所望の積分を数式９の様にM個の積分に分割し各空間Aｊに割り当てる。
【００２１】
【数９】
2014182418.tif 000011
【００２２】
ｇｊは被積分関数をM個の空間へ分配する比であり、数式１０を満たし、特異点ｑｊ，０の近傍でｇｊ（ｑ）〜１、ｇｋ（ｑ）〜０を満たす様な関数である。（ここでｋはｊと異なる。）
【００２３】
【数１０】
2014182418.tif 000012
【００２４】
標本点にたいする重みＷｊ，ｉを数式１１の様に定義し、数式４の数値積分を数式１２の様に一般化する。
【００２５】
【数１１】
2014182418.tif 000013
【００２６】
【数１２】
2014182418.tif 000014
【００２７】
ここでｒｊ，ｉは空間Bにおける特異点ｑｊ，０から標本点ｑｊ，ｉまでの距離である。こうして標本点ｑｊ，ｉとその重みＷｊ，ｉが与えられ単純な関和によって数値積分が実行できる。
【００２８】
分配比ｇｊとしては様々な関数形が考えられるが、M＝３の場合の一例を示す。数式１３のＷ‘ｊ，ｉを用いて、ｇ１（ｑ１，ｉ）を数式１４の様に定義する。ｒｋ，ｉは特異点ｑｋ，０と標本点ｑｊ，ｉの距離であり、ここでｋ≠ｊでもｋ＝ｊでも良い。言い換えればＷ‘ｋ，ｉは、本来空間Aｊから空間Bに対付けられた標本点であるｑｊ，ｉが、仮に空間Aｋから空間Bに対付けられた標本点であったとした時に与えられるべき重みである。
【００２９】
【数１３】
2014182418.tif 000015
【００３０】
【数１４】
2014182418.tif 000016
【００３１】
ｇ２（ｑ１，ｉ）やｇ３（ｑ１，ｉ）は、添え字１，２，３を循環的に置換すると得られる。これらの関数は数式１０を満たす事は明らかであり、また原子１の原子核近傍ではｗ１，ｉが特に小さくなるのでｇ１〜１、ｇ２〜０、ｇ３〜０となりｇｊに求められる性質を満たす事が分かる。なお数式１４の定義から数式１５が導かれる。これを一般化すると標本点ｑｊ，ｉの重みｗｊ，ｉは数式１６で定義される。
【００３２】
【数１５】
2014182418.tif 000017
2014182418.tif 000018
【００３３】
【数１６】
2014182418.tif 000019
【００３４】
以上の一般のMに関する説明は、前述のM=1の場合の自然な拡張になっており、M=1の場合を含む。
【００３５】
分子の分子軌道や全エネルギーについては、以上の方法で得られる標本点やその重みを用いて数値積分を行っても、必要に応じて標本点数を増やせば精度の高い計算結果を得る事が出来る。しかし分子の平衡構造を求める等の目的から、分子軌道法により原子核に働く力（＝エネルギー勾配）を計算する必要が生じる事がある。非特許文献３には数値積分法を用いて原子に働く力（エネルギー勾配）を求め、分子の構造を最適化するための方法が説明されている。ところが非特許文献３によると、従来の方法による標本点やその重みを用いて数値積分を行ったのでは得られた力の精度が不十分な場合が多く、特に周期律表の第４周期以上の原子を含む分子に対しては意味のある結果を得る事が出来なかった。また第２周期以下の原子のみを含む分子の場合でも、多数の標本点を用いる必要があり計算に多大の時間を要した。
【００３６】
原子に働く力の場合には、主要項は注目している原子の原子核に他の原子核や電子から働くクーロン力である。他の原子核から働く力は点電荷間のクーロン力なので正確に計算できる。電子から働くクーロン力は、例えばｚ成分の場合、数式１７の形で表される。
【００３７】
【数１７】
2014182418.tif 000020
【００３８】
被積分関数が数式１７の様な場合、変数変換された数値積分の式（数式４や数式１２）において、ｒ→０でｄｗ/ｄｖ（数式８や数式１１）が実質的にｒ２となるため、数式１７の分母のｒ２と相殺されてしまう。その上通常φ２は原子核近傍で大きくなる。従って数式４や数式１２の右辺の原子核近傍の標本点に対応する項は特に大きくなり支配的な誤差を生じたと考えられる。その他にも（万有）引力の様に強い中心力が働く場合等、数式３のｎが２や２以上に大きな関数を数値積分する場合には、従来の標本点やその重みでは計算が困難、または計算時間が著しく長くなるので改善が求められていた。
【先行技術文献】
【非特許文献】
【００３９】
【非特許文献１】津田孝夫著  「モンテカルロ法とシミュレーション＜三訂版＞」培風館  １９９５年
【非特許文献２】足立裕彦著  「量子材料化学入門」三共出版  １９９１年
【非特許文献３】中川克己  「日本コンピュータ化学会  年報Vol.１１」  ２０１２年  １９４〜２０１頁
【発明の概要】
【発明が解決しようとする課題】
【００４０】
解決しようとする課題は、分子軌道法等の数値シミュレーションにおいて、値（絶対値）が激しく変化する多変数関数が精度良く数値積分できない問題である。
【課題を解決するための手段】
【００４１】
本発明は、仮想的な空間Aにおける標本点から被積分関数が定義された空間Bにおける標本点への変換式を従来の変換式より改善して、各標本点における関数値と標本点の重みの積のバラツキを抑えることを最も主要な特徴とする。
【発明の効果】
【００４２】
各標本点における関数値と標本点の重みの積のバラツキを抑えることにより、関数値が特に大きい特定の標本点における誤差が積分値全体に対して支配的になる事を防げるため、標本点数の増加を積分値の精度向上に直結出来る利点がある。
【図面の簡単な説明】
【００４３】
【図１】図１はGeSe分子の各原子に働く力を両原子の核間距離の関数として示したものである。（実施例２の比較例）
【図２】図２はGeSe分子の各原子に働く力を両原子の核間距離の関数として示したものである。（実施例２）
【図３】図３は本発明の数値積分法を実施するためのプログラムの概要を示したものである。（実施例４）
【発明を実施するための形態】
【００４４】
本発明者は、数値積分に使用する標本点やその重みの与え方を、被積分関数の性質に適応させる必要があると考えた。始めに原子の場合（M=1）を考える。数式１７の様な関数を適切に積分するには、数式８を数式１８の形に書いた時、ｒ→０で１より大きな値に緩慢に収束している右辺のｈ（ｒ）を、ｒ→０の時にφ２の増加を相殺できる程度に１以下に急速に減少させるべきである。
【００４５】
【数１８】
2014182418.tif 000021
【００４６】
ｒ→０でｈ（ｒ）が収束する値ｈ０が小さいほどｈ（ｒ）の減少は早いと考えられる。（ｄｗ/ｄｖは重みの比なのでｈ０は０か又は正の値）さらにｈ０の最適値は注目している原子のZに依存すると予想されるので、ｈ０は調整可能とすべきである。また原子核から必要以上の遠方で標本点が形成されないためには、ｄｗ/ｄｖが数式８の様にｒ→∞で十分に大きくなる事が望ましい。以上の観点より、本発明者は数式８に代わる重みの比の方程式として数式１９に思い至った。
【００４７】
【数１９】
2014182418.tif 000022
【００４８】
ここでR０はｈ０を調整するためのパラメーターである。ｈ（ｒ）は、R０が十分小さな正の値ならｒ→０で急速に減少してｈ０は１より小さな正の値となり、R０＝０ならｈ０＝０となる。この時、数式６で与えられていた標本点の変換の式は数式２０に、数式７で与えられていた係数Γは数式２１に変更する必要がある。以上が本発明の基本的な考え方であり、複数の原子からなる分子等、特異点がM個ある一般の場合も、背景技術に記載した従来の方法に以上の変更点を適合させれば良い。
【００４９】
【数２０】
2014182418.tif 000023
【００５０】
【数２１】
2014182418.tif 000024
【００５１】
即ち、まず各特異点ｑｊ，０に対応する空間Aｊを考え、Aｊにおいて所定の方式に従って標本点ｐｊ，ｉを配置し重みｖｊ，ｉを与える。（ｊ＝１，２，・・・，Ｍ，ｉ＝１，２，・・・，Ｎｊ、Ｎｊは一般にｊに依存する。）その方式としては準乱数を使っても良いが、それに限らない。次に数式２０に従って標本点ｐｊ，ｉから空間Bに変換を行う。各空間においてパラメーターR０や係数Γは異なっていても良い。次に変換で得られた座標を特異点ｑｊ，０の座標分だけ平行移動し標本点ｑｊ，ｉに対応付ける。
【００５２】
さらに数式１０の性質を満たす被積分関数の分配比ｇｊを用いて所望の積分を数式９の様にM個の積分に分割し、各空間Aｊに割り当てる。標本点ｑｊ，ｉの重みｗｊ，ｉは数式２２の様に表される。ここでｒｊ，ｉは空間Bにおける特異点ｑｊ，０から標本点ｑｊ，ｉまでの距離である。
【００５３】
【数２２】
2014182418.tif 000025
【００５４】
ｇｊとしては、数式２３のｗ‘ｋ，ｉを用いて数式１４で与えられる分配比が好適に使用できる。このｗ‘ｋ，ｉを使えば、標本点の重みｗｊ，ｉを数式１６で直接定義できる。ｒｋ，ｉは特異点ｑｋ，０と標本点ｑｊ，ｉの距離であり、ここでｋ≠ｊでもｋ＝ｊでも良い。
【００５５】
【数２３】
2014182418.tif 000026
【００５６】
こうして得られた標本点ｑｊ，ｉとその重みｗｊ，ｉは従来の場合と全く同様に数式１２の数値積分の式の中で使用できる。従って、既存の数値シミュレーションのプログラムにおいても、本発明の方法で得られた標本点と重みを読み込める様にするだけで、精度の高い計算結果が得られる点が本発明の一つの特徴である。
【実施例１】
【００５７】
本実施例においては、孤立したＧｅ原子（３重項状態）に働く力（エネルギー勾配）を計算する。もとより孤立した原子に外部から力が働くはずはないが、逆に計算によって０でない力が得られた場合はその値を誤差そのものと見なす事が出来る。計算の基本となるＧｅ原子の分子軌道は非特許文献２に記載されたDV-Xα法によって求めた。この方法では背景技術に記載された準乱数を用いて空間Aにおける標本点を生成し、各標本点を数式６に従って空間Bへ対応付けて得られた標本点とその重みを用いて分子軌道を求めている。次いでこの分子軌道を用いて非特許文献３に記載された方法で原子に働く力を計算した。こうして得られた結果を比較例とする。
【００５８】
一方空間Aにおける同じ標本点を用いながらも、この標本点を数式２０に従って対応付けて得た空間Bの標本点と数式２２に従って得られた重みを用いてDV-Xα法で分子軌道を求め、次いでこの分子軌道を用いて非特許文献３に記載された方法で原子に働く力を計算し実施例とした。その際、R０を７水準に変化させその効果を見た。すべての計算で１８８００個の標本点を用いている。こうして得られた原子に働く力のベクトルの長さが、（Hartree）/（原子単位の長さ）を単位として表１に示されている。実施例２で示される様に、GeSe分子の核間距離を平衡値より１０％伸ばした時に原子に働く引力がこの単位で０．１程度である。比較例の値が示す様な誤差があっては、分子内の原子に働く力を論ずることが出来ないのは明らかである。本実施例の方法により原子に働く力（＝誤差）は大幅に小さくなった。さらに誤差はR0の値に依存し、本実施例の場合はR0=０．００２が最適と判断された。
【００５９】
【表１】
2014182418.tif 000027
【実施例２】
【００６０】
本実施例においては、主軸をZ方向に向けて配置されたGeSe分子のGe原子（下側）とSe原子（上側）に働く力を、両原子間の核間距離をパラメーターとして計算した。準乱数を用いて空間A１においてGe原子用の標本点を生成し、同様に空間A２においてSe原子用の標本点を生成した。各標本点を数式６に従って空間Bへ対応付けて得た７５２００個の標本点とその重みを用いて、DV-Xα法で分子軌道を求めた。この分子軌道を用いて非特許文献３に記載された方法に従ってエネルギー勾配を計算した。こうして得られた結果を比較例として図１に示す。ここで横軸１０１は実測の平衡核間距離（＝２．１３４６Å）を１．０とした両原子の核間距離であり、縦軸１０２は実施例１と同じ単位で表した両原子に働く力のｘ方向、y方向、ｚ方向の成分である。
【００６１】
空間A１と空間A２において比較例と同様にして標本点を生成し、これからR０=０．００２として数式２０に従って空間へ対応付けて得た７５２００個の標本点を用いて、DV-Xα法で分子軌道を求めた。この分子軌道を用いて非特許文献３に記載された方法に従ってエネルギー勾配を計算した。こうして得られた結果を実施例として図１と同じ要領で図２に示す。
【００６２】
図２においては、核間距離２０１が短いと下側に配置されたGe原子には下向きの力（２０７）が働き、上側に配置されたSe原子には絶対値が同じ上向きの力（２０８）が働き（斥力）、両原子間の距離が長いと下側に配置されたGe原子には上向きの力（２０７）が働き、上側に配置されたSe原子には絶対値が同じ下向きの力（２０８）が働く（引力）事が分かり、合理的な結果となっている。両原子に働く力がともに０となるのは核間距離＝０．９６で、実測の平衡距離１．００よりやや短いが、同時に計算された全エネルギー核間距離＝０．９６で極小となる事から、実測値との不一致は数値積分法の誤差によるものではないと推測される。一方図１においては、大よその傾向は図２と似ているものの、各核間距離で両原子に働くｚ方向の力の絶対値は異なり、しかも本来力が働くはずがないｘ方向やｙ方向にも大きな力が働くと言う不合理な結果となった。本発明の数値積分法による精度向上の効果は明らかである。
【実施例３】
【００６３】
本実施例においては、これまでに説明してきたエネルギー勾配を用いて、正四面体対称のSiCl4分子の構造最適化を行った。使用した構造最適化の方法については非特許文献３に記載がある。まずSiCl４分子に表２の欄１）に示すような意図的に正四面体対称からはずした初期構造を与えた。次に準乱数を用いて空間A１においてSi原子用の標本点を生成し、同様にして空間A２〜空間A５において、４個のCl原子用の標本点を生成した。各標本点を数式６に従って空間Bへ対応付けて得た９４０００個の標本点とその重みを用いて、DV-Xα法で分子軌道を求めた。この分子軌道を用いて非特許文献３に記載された方法に従って各原子に働く力を計算した。各原子に働く力に合わせて各原子を少し移動させ、最初の工程と同様にして標本点と重みを得て分子軌道を求め各原子に働く力を再計算した。この様な工程を１２回繰り返した所で、各原子に働く力に合わせて原子を動かしても全エネルギーが下がらなくなったので計算を終了した。最終の構造を表２の欄２）に示す。
【００６４】
SiCｌ４分子に対して比較例と同じ初期構造を与えた。空間A１〜空間A５において比較例と同様にして標本点を生成し、これからR０=０．００２として数式１１に従って空間Bへ対応付けて得た９４０００個の標本点を用いて、DV-Xα法で分子軌道を求めた。この分子軌道を用いて非特許文献３に記載された方法に従って各原子に働く力を計算した。以下、比較例と同様にして構造最適化の工程を１１回繰り返した所で、各原子に働く力に合わせて原子を動かしても全エネルギーが下がらなくなったので計算を終了した。最終の構造を表２の欄３）に示す。あわせて実測によるSiCi４分子の構造を表２の欄４）に示す。
【００６５】
【表２】
2014182418.tif 000028
【００６６】
実施例の最終構造では、Si-Clの長さは実測値に誤差１％以内で一致し、∠ClSiClは誤差±２°以内で一致した。一方比較例の最終構造では、Si-Clの長さの誤差は＋８％に及び、∠ClSiClの誤差は−３０°〜＋１０°に及んだ。比較例ではエネルギー勾配の計算精度が不十分なため、構造最適化の工程において適切な方向に構造を変化させる事が出来なかったためと思われる。本発明の方法による積分精度向上の効果は明らかである。
【実施例４】
【００６７】
本発明の方法で標本点とその重みを生成し数値積分を行うプログラムをFortran言語で作成した。その概要を図３に示す。ここでMは分子を構成する原子の数である。Nｊは原子ｊに割り振る標本点の数であり、ｊに異存していてもよい。プログラムの３０１〜３１７で標本点とその重みを求める。計算の具体的な内容についてはこれまで説明してきた数式を引用して簡略に記述している。本実施例では被積分関数の分配比ｇｉとして数式１４の関数を用いている。３１８〜３２４では得られた標本点とその重みを用いて数値積分を実行している。数値積分は全体を通して実行されるが、本発明による改良の本質的な部分は３０１〜３１７にある。しかし数値積分法として一体となったプログラムとする場合は、３２１を３１５の後に挿入しても良い。図５から明らかなように、本実施例のプログラムは単純で計算の負荷は少ないと思われるが、大きなモデルの計算を行う等の理由で計算時間の短縮が求められる場合は、例えば３０６〜３１６のｄｏループをOpenMP等の手法により並列化して複数のコアを有する計算機にかければよい。
【００６８】
実施例１から実施例３におけるエネルギー勾配の計算においては、５０１〜５１７行のプログラムで得た標本点とその重みをDV−Xα法によるプログラムに読み込ませて所望の原子または分子の分子軌道を求めた。続いて非特許文献３に記載された方法でエネルギー勾配を求める既存のプログラムを全く改変なく使用して原子に働く力を求めた。既存の数値計算法を含むプログラムに本発明のプログラムの結果を渡すだけで大幅な計算精度の向上が図れる点も本発明の特徴である。
【００６９】
以上の実施例を通して、本発明の標本点とその重みを用いれば、値が大幅に変化する多変数の関数を、従来の方法よりはるかに精度良く積分出来る事を示した。一方従来の方法でも積分が可能だった関数も、標本点の数を少なくしても同等の精度で積分できるので、計算時間を短縮する事も可能となる。
【産業上の利用可能性】
【００７０】
クーロン力や万有引力等、空間で値が急激に変化する量の積分を伴う数値シミュレーションにおいて、これまで数値積分の誤差のため扱いが困難だったモデルが扱える、及び/又は必要な計算時間が大幅に短縮出来ると期待される。
【符号の説明】
【００７１】
１０１、２０１    両原子の核間距離
１０２、２０２    各原子に働く力
１０３、２０３    Ge原子に働く力のｘ成分
１０４、２０４    Se原子に働く力のｘ成分
１０５、２０５    Ge原子に働く力のｙ成分
１０６、２０６    Se原子に働く力のｙ成分
１０７、２０７    Ge原子に働く力のｚ成分
１０８、２０８    Se原子に働く力のｚ成分
３０１〜３１７    本発明の標本点とその重みの生成方法を実施するためのプログラムのステートメント
３１８〜３２４    本発明により与えられた標本点とその重みを用いて数値積分を行うプログラムのステートメント
【図１】
2014182418.tif 000033
【図２】
2014182418.tif 000034
【図３】
2014182418.tif 000035