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