特許 5852032 多変数関数の数値積分用標本点とその重みの生成方法及び生成プログラム

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

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

(11)【特許番号】5852032

(24)【登録日】2015年12月11日

(45)【発行日】2016年2月3日

(54)【発明の名称】多変数関数の数値積分用標本点とその重みの生成方法及び生成プログラム

(51)【国際特許分類】

   G06F 17/11 20060101AFI20160114BHJP

   G06F 19/00 20110101ALI20160114BHJP

【ＦＩ】

   G06F17/11

   G06F19/00 110

【請求項の数】6

【全頁数】21

(21)【出願番号】特願2013-54602(P2013-54602)

(22)【出願日】2013年3月18日

(65)【公開番号】特開2014-182419(P2014-182419A)

(43)【公開日】2014年9月29日

【審査請求日】2014年9月30日

(73)【特許権者】

【識別番号】312010799

【氏名又は名称】中川 克己

(72)【発明者】

【氏名】中川 克己

【審査官】 田中 幸雄

(56)【参考文献】

【文献】 特開２００７−１０９０２３（ＪＰ，Ａ）

【文献】 特開２００８−０２１２５９（ＪＰ，Ａ）

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

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

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

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

Ｇ０６Ｆ １７／１１

Ｇ０６Ｆ １９／００

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

【請求項1】

コンピューターに
空間Bで定義された被積分関数の特異点ｑ_ｊ，０（ｊ＝１，２，・・・，Ｍ）を取得させる工程と、
特異点ｑ_ｊ，０毎に対応する空間A_ｊを設定させ、
各空間A_ｊにおいて、
所定の方式でN_ｊ個の標本点ｐ_ｊ，ｉ（ｉ＝１，２，・・・，Ｎ_ｊ）を配列させる工程と、
標本点ｐ_ｊ，ｉに対して重みｖ_ｊ，ｉを与えさせる工程と、
空間A_ｊから空間Bへの標本点の変換式を設定させる工程と、
標本点ｐ_ｊ，ｉを該変換式により変換しさらに空間Bの原点と特異点ｑ_ｊ，０の位置の差分だけ移動して空間Bの標本点ｑ_ｊ，ｉに対応付けさせる工程と、
空間Bにおける標本点ｑ_ｊ，ｉの重みと空間A_jにおける標本点ｐ_ｊ，ｉの重みの比ｄｗ/ｄｖを該変換式に基づいて算出し、
被積分関数の空間A_ｊへの分配比ｇ_ｊを標本点ｑ_ｊ，ｉと各特異点ｑ_k，０ (k＝１，２，・・・，Ｍ）の位置関係に基づいて算出し、
標本点ｑ_ｊ，ｉの重みＷ_ｊ，ｉとして

【数25】

を算出させる工程と、
標本点ｑ_ｊ，ｉの位置とその重みＷ_ｊ，ｉをメモリに保存または出力させる工程とを実行させる多変数関数の数値積分用の標本点とその重みの生成方法において、
空間A_ｊの標本点ｐ_ｊ，ｉは面心立方格子状に配列される事を特徴とする多変数関数の数値積分用の標本点とその重みの生成方法。

【請求項2】

請求項１に記載の多変数関数の数値積分用の標本点とその重みの生成方法において、
前記重みの比ｄｗ/ｄｖは、特異点ｑ_ｊ，０と標本点ｑ_ｊ，ｉとの距離ｒに対してｒ^２ｈ（ｒ）の形で依存し、
関数ｈ（ｒ）はｒ→０に対して減少しつつ所定の値に収束する関数であり、
該所定の値はｈ（ｒ）が含むパラメーターの調整によって０または任意の正の値に出来るものである事を特徴とする、
多変数関数の数値積分用の標本点とその重みの生成方法。

【請求項3】

請求項１および請求項２に記載の多変数関数の数値積分用の標本点とその重みの生成方法において、
前記空間A_ｊの標本点ｐ_ｊ，ｉから対応付けられた標本点ｑ_ｊ，ｉは、前記空間Bにおいて空間A_ｊに対応する特異点ｑ_ｊ，０を含む少なくとも一つの平面に対して面対称に配列される事を特徴とする多変数関数の数値積分用の標本点とその重みの生成方法。

【請求項4】

コンピューターに
空間Bで定義された被積分関数のM個の特異点ｑ_ｊ，０（ｊ＝１，２，・・・，Ｍ）の座標を入力させ、
各特異点ｑ_ｊ，０に対応する空間A_ｊを設定させ、
空間A_ｊにおいて、
空間A_ｊから空間Bへの標本点の変換式を設定させ、
空間A_ｊにおける標本点数を設定させ、
所定の方式に従って標本点ｐ_ｊ，ｉ（ｉ＝１，２，・・・，Ｎ_ｊ）を生成させ、
標本点ｐ_ｊ，ｉに付随する重みｖ_ｊ，ｉを設定させ、
前記変換式に従って標本点ｐ_ｊ，ｉを変換しさらに空間Bの原点と特異点ｑ_ｊ，０の位置の差分だけ移動し空間Bの標本点ｑ_ｊ，ｉに対応付けさせ、
空間Bにおける標本点ｑ_ｊ，ｉの重みと空間A_jにおける標本点ｐ_ｊ，ｉの重みの比ｄｗ/ｄｖを該変換式に基づいて算出させ、
被積分関数の空間A_ｊへの分配比ｇ_ｊを標本点ｑ_ｊ，ｉと各特異点ｑ_k，０ (k＝１，２，・・・，Ｍ）の位置関係に基づいて算出させ、
標本点ｑ_ｊ，ｉの重みＷ_ｊ，ｉとして数式２６を算出し、

【数26】

標本点ｑ_ｊ，ｉとその重みＷ_ｊ，ｉをメモリに保存または出力させる多変数関数の数値積分用の標本点とその重みの生成プログラムにおいて、
空間A_ｊの標本点ｐ_ｊ，ｉは面心立方格子状に配列される事を特徴とする多変数関数の数値積分用の標本点とその重みの生成プログラム。

【請求項5】

請求項４に記載の多変数関数の数値積分用の標本点とその重みの生成プログラムにおいて、
前記重みの比ｄｗ/ｄｖは、空間Bの特異点ｑ_ｊ，０と標本点ｑ_ｊ，ｉとの距離ｒに対してｒ^２ｈ（ｒ）の形で依存し、
関数ｈ（ｒ）はｒ→０に対して減少しつつ所定の値に収束する関数であり、
該所定の値はｈ（ｒ）が含むパラメーターの調整によって０または任意の正の値に出来るものである事を特徴とする、
多変数関数の数値積分用の標本点とその重みの生成プログラム。

【請求項6】

請求項４および請求項５に記載の多変数関数の数値積分用の標本点とその重みの生成プログラムにおいて、
前記空間A_ｊの標本点ｐ_ｊ，ｉから対応付けられた標本点ｑ_ｊ，ｉは、前記空間Bにおいて空間A_ｊに対応する特異点ｑ_ｊ，０を含む少なくとも一つの平面に対して面対称に配列される事を特徴とする多変数関数の数値積分用の標本点とその重みの生成プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、数値シミュレーション等において必要とされる多変数関数の数値積分を高精度で行うために好適な標本点と重みの生成方法及びその生成のためのプログラムに関するものである。

【背景技術】

【0002】

数値シミュレーション等において、必要とされる関数の積分が解析的には不可能または困難な場合がある。その場合でも被積分関数が１変数の関数であれば、Ｓｙｍｐｓｏｎの公式を始めとする各種公式を用いて任意の精度で数値積分を行う事が出来る。しかし被積分関数が多変数関数の場合は、単純に重積分を行ったのでは計算量が著しく増大するので、モンテカルロ法が用いられる事が多い。モンテカルロ法とは、非特許文献１に記載されている様に空間の点ｐにおける被積分関数をf（ｐ）とする時、数式１に示す積分を数式２に示す形で近似する方法である。

【0003】

【数1】

【0004】

【数2】

【0005】

ここでp_iは積分を行う空間に分布するN個の標本点の一つで、例えばｆ（p）が３変数の関数ならｐ_ｉ＝（ｓ_ｉ，ｔ_ｉ，ｕ_ｉ）の様に３個の座標値を持つ。またv_ｉは標本点p_iに割り当てられる重み（空間要素の体積）である。標本点としては準乱数の組が広く用いられている。準乱数とは、適当な無理数を整数倍した値の小数点以下の部分であり、整数として１〜Nを順次与える事によってN個の準乱数が得られる。また３変数の関数の場合は、例えば√２、√３、√５の様な相異なる３個の無理数に対し整数として１〜Nを順次与える事によって標本点がN個得られる。この様にして得られた標本点は空間の中で一様に近い分布をする。通常は完全に一様な分布と見なしてｖ_ｉ＝Ｖ／Ｎを与える。ここでVは積分を行う空間の体積である。但し積分を行う空間の表面に位置する標本点には、これと異なる値を与えても良い。数式２の誤差は概ねNに反比例するので、高い精度が必要なら通常多数の標本点を使う必要がある。

【0006】

数値積分が必要となる関数の中には、積分を行う空間の中でその値（絶対値）を大きく変えるものがある。例えば数値積分が必要となる典型的な数値シミュレーションである分子軌道法の場合、被積分関数として数式３の形の関数が現れる事が多い。

【0007】

【数3】

【0008】

ここでＺは分子を構成する原子の原子番号、ｒは原子核と標本点との距離、ｎは通常１か１以上である。例えば分子の全エネルギーを求めるには原子核と電子のクーロン相互作用のエネルギーを計算する必要があるが、その場合はｎ＝１である。１/ｒ^ｎの寄与に加えて、波動関数φ_１及びφ_２は原子核に近づくにつれ値が大きくなるので数式３は原子核近傍で著しく大きな値を取る事が多い。原子核の位置の様に被積分関数の値が著しく大きくなる点を以下では被積分関数の特異点、又は単に特異点と呼ぶ。一方原子核から遠くなると一般に波動関数は０に向かって指数関数的に収束する。そのため数式２の右辺の各項の大きさは著しく不揃いとなり、大きな値をとる項の誤差が支配的となってNを大きくしてもその割には積分の精度が向上しない。この効果は特にZの大きな原子を含む分子の場合に顕著となる。

【0009】

この問題を改善するため変数の変換が利用される事が多い。始めに原子が1個だけある場合について説明する。始めに仮想的な空間Aを考え、この空間における標本点p_iに対し、適当な変換式によって被積分関数が定義された空間Bにおける標本点ｑ_ｉを対応させる。所望の積分は空間Bにおいて数式４の形の数値積分で近似される。

【0010】

【数4】

【0011】

ｗ_iは標本点ｑ_ｉの重みである。空間Aと空間Bの標本点の重みの比ｄｗ/ｄｖを使うとｗ_iは数式５の様に表される。

【0012】

【数5】

【0013】

ここで変換式を調整して、空間Bの標本点を、関数ｆ（ｑ）の値が大きい場所では密に分布させ、ｆ（ｑ）の値が小さい場所では疎に分布させることにする。するとｆ（ｐ）の値が大きい場所ではｄｗ/ｄｖは小さくなり、ｆ（ｐ）の値が小さい場所ではｄｗ/ｄｖが大きくなるため、数式４の右辺の各項の値のバラツキが小さくなる。このため数式４を直接数値積分するより、Nが同じでも精度が高くなる。見方を変えると、空間Bにおける被積分関数ｆ（ｐ）を、空間Aにおいて値の変化が少ない関数ｆ（ｐ）ｄｗ/ｄｖに変換して数値積分すると考える事も出来る。

【0014】

３変数の場合について具体的な例を考えてみる。非特許文献２にならって空間Aとしてｓ＝[０，１]、ｔ＝[０，１]、ｕ＝[０，１]の立方体を想定し、この中に準乱数を使って標本点を分布させる。空間Bの座標系として極座標系を採用し、数式１５の変換式に従って標本点p_i（ｓ_ｉ，ｔ_ｉ，ｕ_ｉ）を標本点ｑ_i（Φ_ｉ，θ_ｉ，ｒ_ｉ）に対応付ける。

【0015】

【数6】

【0016】

ここでR_０は標本点の分布の広がりを調整するパラメーターである。Γはｕ→１とする時ｒ→∞となるように定められる係数であり数式７で与えられる。この時ｄｗ/ｄｖは数式８の様に与えられる。

【0017】

【数7】

【0018】

【数8】

【0019】

従って関数の特異点（原子核の位置）の近傍で被積分関数ｆ（ｑ）が数式３に従って大きくなっても、ｎ＝１であれば数式４の右辺各項は特に大きくならないので数値積分の精度は低下し難い。

【0020】

以下ではM個の原子を含む一般の分子について説明する。この場合、空間BにM個の特異点ｑ_ｊ，０（原子核の位置，ｊ＝１，２，・・・，Ｍ）を含む事になる。空間Aの標本点からM個の特異点を含む空間Bの標本点に、単一の変換式で対応づけられると理想的であるが、今のところ適当な変換式が見つかっていない。そこで各特異点に対応する空間A_１、空間A_２、・・・、空間A_Mを考え、各空間において所定の方式に従って標本点ｐ_ｊ，ｉ（ｊ＝１，２，・・・，Ｍ，ｉ＝１，２，・・・，Ｎ_ｊ）を配置することにする。Ｎ_ｊは特異点ｑ_ｊ，０の性質によって最適な値が異なるので一般にはｊに依存する。次に数式６に従って標本点ｐ_ｊ，ｉから空間Bに変換する。各空間によってR_０や係数Γも異なっていて良い。さらに特異点ｑ_ｊ，０の座標分だけ平行移動し標本点ｑ_ｊ，ｉとする。
ここで所望の積分を数式９の様にM個の積分に分割し各空間A_ｊに割り当てる。

【0021】

【数9】

【0022】

ｇ_ｊは被積分関数をM個の空間へ分配する比であり、数式１０を満たし、特異点ｑ_ｊ，０の近傍でｇ_ｊ（ｑ）〜１、ｇ_ｋ（ｑ）〜０を満たす様な関数である。（ここでｋはｊと異なる。）

【0023】

【数10】

【0024】

標本点にたいする重みＷ_ｊ，ｉを数式１１の様に定義し、数式４の数値積分を数式１２の様に一般化する。

【0025】

【数11】

【0026】

【数12】

【0027】

ここでｒ_ｊ，ｉは空間Bにおける特異点ｑ_ｊ，０から標本点ｑ_ｊ，ｉまでの距離である。こうして標本点ｑ_ｊ，ｉとその重みＷ_ｊ，ｉが与えられ単純な積算によって数値積分が実行できる。

【0028】

分配比ｇ_ｊとしては様々な関数形が考えられるが、M＝３の場合の一例を示す。数式１３のＷ‘_ｋ，ｉを用いて、ｇ_１（ｑ_１，ｉ）を数式１４の様に定義する。ｒ_ｋ，ｉは特異点ｑ_ｋ，０と標本点ｑ_ｊ，ｉの距離であり、ここでｋ≠ｊでもｋ＝ｊでも良い。言い換えればＷ‘_ｋ，ｉは、標本点ｑ_ｊ，ｉが、仮に空間A_ｋから空間Bに対付けられた標本点であったとした時に与えられるべき重みである。

【0029】

【数13】

【0030】

【数14】

【0031】

ｇ_２（ｑ_１，ｉ）やｇ_３（ｑ_１，ｉ）は、添え字１，２，３を循環的に置換すると得られる。これらの関数は数式１０を満たす事は明らかであり、また原子１の原子核近傍ではｗ_１，ｉが特に小さくなるのでｇ_１〜１、ｇ_２〜０、ｇ_３〜０となりｇ_ｊに求められる性質を満たす事が分かる。なお数式１４の定義から数式１５が導かれる。これを一般化すると標本点ｑ_ｊ，ｉの重みｗ_ｊ，ｉは数式１６で定義される。

【0032】

【数15】

【0033】

【数16】

【0034】

以上の一般のMに関する説明は、前述のM=1の場合の自然な拡張になっており、M=1の場合を含む。

【0035】

分子の分子軌道や全エネルギーについては、以上の方法で得られる標本点やその重みを用いて数値積分を行っても、必要に応じて標本点数を増やせば精度の高い計算結果を得る事が出来る。しかし分子の平衡構造を求める等の目的から、分子軌道法により原子核に働く力（＝エネルギー勾配）を計算する必要が生じる事がある。非特許文献３には数値積分法を用いて原子に働く力（エネルギー勾配）を求め、分子の構造を最適化するための方法が説明されている。ところが非特許文献３によると、従来の方法による標本点やその重みを用いて数値積分を行ったのでは得られた力の精度が不十分な場合が多く、特に周期律表の第４周期以上の原子を含む分子に対しては意味のある結果を得る事が出来なかった。また第２周期以下の原子のみを含む分子の場合でも、多数の標本点を用いる必要があり計算に多大の時間を要した。

【0036】

原子に働く力の場合には、主要項は注目している原子の原子核に他の原子核や電子から働くクーロン力である。他の原子核から働く力は点電荷間のクーロン力なので正確に計算できる。電子から働くクーロン力は、例えばｚ成分の場合、数式１７の形で表される。

【0037】

【数17】

【0038】

被積分関数が数式１７の様な場合、特異点の近傍で関数の値は正負を挟んで激しく変化する。標本点の空間要素内でも関数の変化が大きく標本点での値が空間要素内での値を適切に代表しなくなり、モンテカルロ法の様な従来の方法では積分で高い精度を得るのが困難、または計算時間が著しく長くなるので改善が求められていた。

【先行技術文献】

【非特許文献】

【0039】

【非特許文献1】津田孝夫著 「モンテカルロ法とシミュレーション＜三訂版＞」培風館 １９９５年

【非特許文献2】足立裕彦著 「量子材料化学入門」三共出版 １９９１年

【非特許文献3】中川克己 「日本コンピュータ化学会 年報Vol.１１」 ２０１２年 １９４〜２０１頁

【非特許文献4】中川克己 「DV-Xα研究協会会報 Vol.２３」 ２０１０年 ８８〜９３頁

【発明の概要】

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

【0040】

解決しようとする課題は、分子軌道法等の数値シミュレーションにおいて、値（絶対値）が特に激しく変化する多変数関数が精度良く数値積分できない問題である。

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

【0041】

本発明は、被積分関数が定義された空間Bにおける標本点を、空間Aにおいて面心立方格子状に規則的に配列された標本点から、変数変換によって空間Bに対応付けて生成する事を最も主要な特徴とする。

【発明の効果】

【0042】

空間Aの標本点を規則的に配列する事により、標本点における被積分関数の値を内挿する２次式が得られる。得られた２次式を積分すれば、変化の激しい関数でも高い精度で積分出来る利点がある。

【図面の簡単な説明】

【0043】

【図1】図１は面心立方格子の１つの格子点とその１２個の最近接格子点、及び中心の格子点の空間要素を示したものである。

【図2】図２は準乱数による従来の標本点と、面心立方格子をなす標本点とを用いてモデル関数の数値積分を行った際に生じた誤差の比較するものである。

【図3】図３はSnTe分子の各原子に働く力を両原子の核間距離の関数として示したものである。（実施例２の比較例）

【図4】図４はSnTe分子の各原子に働く力を両原子の核間距離の関数として示したものである。（実施例２）

【図5】図５は本発明の数値積分法を実施するためのプログラムの概要を示したものである。（実施例４）

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

【0044】

本発明者は、正確な積分値を得るためには標本点の重みを正確に与える必要があると考えた。重みｗ_ｊ，ｉは数式１１に示される３つの成分から構成される。この内、ｇ_ｊやｄｗ/ｄｖは解析的な形で与えられるので誤差を生じ難いが、ｖ_ｊ，ｉには問題がある。空間Aにおける標本点ｐ_ｊ，ｉを生成した際、通常その重みとしてV/N（＝空間の体積/標本点数）を与えるが、これは標本点がVの体積の中に完全に均一に分布している限り正しい。しかし準乱数を用いて標本点を生成すると、標本点の生成に使用した無理数の組にもよるが、実際の標本点の分布にはかなりの不均一が生じている。従って空間Aにおける標本点の分布を完全に均一にする事によって数値積分の精度を向上し得るはずである。

【0045】

本発明者は標本点の完全に均一な分布の例として、Au、Ag、Cu等の結晶格子として知られる面心立方格子（FCC）を取り上げた。図１に示す様に、FCCでは任意に選択した格子点（中心格子点１０１）の周囲には１２個の最近接格子点１０２がある。（図１では４個の最近接格子点のみが矢印で指示されているが、白丸で示す１２個の格子点は全て最近接格子点である。）非特許文献４には、空間における任意の関数ｆを３変数の２次式で内挿し、中心格子点の近傍での２次式の関数形を中心格子点における関数値ｆ（ｐ_０）と１２個の最近接格子点における関数値ｆ（ｐ_i）で決定する方法が記載され、この２次式を解析的に積分する事により、中心格子点１０１が占める空間要素１０３における関数ｆの積分値として数式１８が与えられている。ここでｄはFCCの格子定数の１/２の長さである。また空間要素の体積は２ｄ^３である。また積分を行う体積が十分に大きければ２ｄ^３＝V/Nとして良い。

【0046】

【数18】

【0047】

即ち１２個の最近接格子点がある場合、ｆ（ｐ_０）は自分自身の空間要素での積分で１回、最近接格子点の空間要素の積分で１２回使われるので、積分全体への寄与は合計ｆ（ｐ_０）２ｄ^３となり、中心格子点の重みとして結局V/Nを用いて良い事になる。これは２次式による内挿をした上で得られた結果なので、標本点を正確にFCCの格子点の様に配列すれば、V/Nが２次の精度を与える重みとなる事を意味し、特に変化の激しい関数を数値積分する際には効果を発揮する。

【0048】

その効果を確認するため、準乱数を用いた標本点を用いる従来の方法と、FCCの格子点を標本点とする方法の二通りの方法で数式１９の数値積分を行った。

【0049】

【数19】

【0050】

ここで分母は規格化のための係数で、十分に広い空間で積分を行えばI＝１．０となる。実際にはａ＝０．１としてｘ＝[−１、１]、ｙ＝[−１、１]、ｚ＝[−１、１]の範囲で、標本点の数N＝４３２からN＝２１２９６まで９水準で積分を行った。図２に両者が示した誤差２０２（＝ I -１．０）を示す。標本点数２０１の増加とともに、いずれの誤差も概ね減少していくが、特にFCCの標本点を用いた場合２０４は、準乱数の標本点を用いた場合２０３より急速に誤差が減少していく。N=3456以上ではほぼ一定となるが、これは積分の範囲が有限なため生じた誤差であり数値積分の誤差ではない。

【0051】

数式１２は空間Bにおける関数ｆの数値積分であると同時に、空間Aにおける関数ｆｇ_ｊｄｗ/ｄｖの数値積分と見る事もできる。従って空間Aの標本点の重みの精度を高めると、空間Bにおける被積分関数ｆの数値積分の精度を高める事になる。従って空間AでFCCの標本点を使用すれば空間Bにおける所望の数値積分の精度を高める事ができる。

【0052】

しかし一つ注意を要する点がある。数式１８から分かる様に、積分領域の周辺部に位置する標本点に対しては重み＝V/Nは正確ではない。幸い図２の被積分関数の場合には、周辺部で被積分関数の値が殆ど０であり、周辺部の標本点の重みの不正確さにより生じる誤差は無視しうるのでFCCによる計算例では極めて高い精度が得られた。逆に言えば、周辺部で被積分関数の値が十分に小さくならないとFCCを用いても精度が得られない。そのため背景技術で記載した従来の方法においても、空間Aの周辺部で被積分関数ｆｇ_ｊｄｗ/ｄｖの値が小さくなる様に、重みの比ｄｗ/ｄｖとしてｒ→０で０に収束する数式８の関数を採用したのである。だが被積分関数が数式１７の形をとると、従来のｄｗ/ｄｖではその効果が十分ではなくなる。そこで本発明者は、数式８の右辺の[ ]内の式をh（r）とする時、現状ではｒ→０で１より大きな値に緩慢に収束しているｈ（ｒ）を、０または１よりも小さな所定の正の値に急速に収束する関数に改良した。この様な重みの比ｄｗ/ｄｖの一例として数式２０を示す。

【0053】

【数20】

【0054】

ここでR_０はｄｗ/ｄｖのｒ→０での収束の速さを調整するためのパラメーターである。ｒ→０で、R_０が十分小さな正の値ならｄｗ/ｄｖは急速に減少して１より小さな正の値となり、R_０＝０なら０となる。この時、数式６で与えられていた標本点の変換の式は数式２１に、数式７で与えられていた係数Γは数式２２に変更する。以下、複数の原子からなる分子等、特異点がM個ある一般の場合について説明する。

【0055】

【数21】

【0056】

【数22】

【0057】

即ち、各特異点ｑ_ｊ，０に対応する空間A_ｊにおいて、標本点ｐ_ｊ，ｉをFCCの格子点を構成するように配置し、重みｖ_ｊ，ｉ＝２ｄ^３を与える。（Ｍ，ｉ＝１，２，・・・，Ｎ_ｊ）標本点の数Ｎ_ｊは空間毎に異なっていても良い。次に数式２１に従って標本点ｐ_ｊ，ｉから空間Bに変換を行う。空間A_ｊによってパラメーターR_０や係数Γは異なっていても良い。次に標本点の座標を特異点ｑ_ｊ，０の座標分だけ平行移動し標本点ｑ_ｊ，ｉに対応付ける。

【0058】

さらに数式１０の性質を満たす被積分関数の分配比ｇ_ｊを用いて所望の積分を数式９の様にM個の積分に分割し、各空間A_ｊに割り当てる。標本点ｑ_ｊ，ｉの重みｗ_ｊ，ｉは数式２３の様に表される。

【0059】

【数23】

【0060】

ｇ_ｊとしては、数式２４のｗ‘_ｋ，ｉを用いて数式１４で与えられる分配比が好適に使用できる。ｒ_ｋ，ｉは空間Bにおける特異点ｑ_ｋ，０から標本点ｑ_ｊ，ｉまでの距離である。ここでｋ≠ｊでもｋ＝ｊでも良い。このｗ‘_ｋ，ｉを使えば、標本点の重みｗ_ｊ，ｉを数式１６で直接定義できる。

【0061】

【数24】

【0062】

被積分関数の性質によっては、空間Aにおいて標本点が規則正しく配列される事を利用してさらに積分の精度を高める事が出来る。数式１７において、原子核（特異点）近傍においてFｚの被積分関数は著しく大きな値を持つが、φ^２で大きな寄与をするのはs軌道の成分であるためφ^２はほぼ球対称の性質を持つ。そこにｃｏｓ（θ）がかかるので被積分関数は空間Bにおいて全体としてθ＝９０°の面（ｘｙ平面）に対して奇となって積分値は本来０に近いはずであるが、実際には各標本点において生じる誤差が積み重なって大きな誤差となる。もし空間Bにおいて標本点がθ＝９０°の面（ｘｙ面）に対して面対称になる様に配列されていると、標本点ｑにおいて生じる誤差と、ｑと面対称の位置にある標本点ｑ^＊で生じる誤差とは絶対値が同じで符号が異なり、相殺し合って積分値の誤差としては表れ難くなると期待される。同様に標本点がｙｚ面（φ＝９０°と２７０°の面）に対して面対称に配列されているとFｘの精度が、標本点がｘｚ面（φ＝０°と１８０°の面）に対して面対称に配列されているとFｙの精度が高まると予想される。

【0063】

空間Aにおける標本点の配列を調整する事によって、空間Bにおいて標本点を特異点ｑ_ｊ，０を含む面に対して面対称に配列する事が出来る。空間A_ｊでφに対応付けられる変数ｓの方向の格子点の数を、例えば８にすると、空間Bでφ＝０°、４５°、・・・、３１５°の方向に標本点が配列され、yz面に対しても、ｘｚ面に対しても面対称と出来る。一般に変数ｕの方向の格子点の数を偶数にすればよい。また空間Aでθに対応付けられる変数ｔ＝[０、１]で、ｔ＝０．５に対して対称となる様に標本点を配列すると、空間Bでｘｙ面に対して対称に標本点を配列できる。

【0064】

こうして得られた標本点ｑ_ｊ，ｉとその重みｗ_ｊ，ｉはモンテカルロ法等従来の場合と全く同様に数式１２の数値積分の式の中で使用できる。従って、既存の数値シミュレーションのプログラムにおいても、本発明の方法で得られた標本点と重みを読み込める様にするだけで、数値積分の精度が向上し正確なシミュレーションが可能となる点が本発明の一つの特徴である。

【実施例1】

【0065】

本実施例においては、孤立したSn原子（３重項状態）に働く力（エネルギー勾配）を計算する。もとより孤立した原子に外部から力が働くはずはないが、逆に計算によって０でない力が得られた場合はその値を誤差そのものと見なす事が出来る。計算の基本となるSn原子の分子軌道は非特許文献２に記載されたDV-Xα法によって求めた。その際に３種類の標本点とその重みを用意しDV-Xα法のプログラムに読み込ませた。
比較例１）空間Aの標本点は準乱数を使って生成し、空間Aの標本点と空間Bの標本点の重みの比ｄｗ/ｄｖは数式８を使用した。
比較例２）空間Aの標本点は準乱数を使って生成し、空間Aの標本点と空間Bの標本点の重みの比ｄｗ/ｄｖは数式２０を使用した。
実施例）空間Aの標本点はFCCを使って生成し、空間Aの標本点と空間Bの標本点の重みの比ｄｗ/ｄｖは数式８を使用した。
次いで３種の分子軌道を用いて非特許文献３に記載された方法で原子に働く力を計算した。すべての計算で３７６００個の標本点を用いた。結果を表１に示す。

【0066】

こうして得られた原子に働く力のベクトルの長さが、（Hartree）/（原子単位の長さ）を単位として表１に示されている。実施例２で示される様に、SnTe分子の核間距離を平衡値より２０％伸ばした時に原子に働く引力がこの単位で０．０７程度である。比較例１の値が示す様な誤差があっては、分子内の原子に働く力を論ずることが出来ないのは明らかである。比較例２にはかなりの改善が見られるが満足できる水準ではない。実施例に示す様にFCCを用いる事によって十分な精度が得られた。なお、実施例では空間Bにおいて標本点をｘｙ平面、ｙｚ平面、ｚｘ平面のいずれに対しても面対称となる様に配置した。実施例で非常に精度が高い結果が得られたのにはその寄与もある。標本点の面対称配列がこれほどの効果を持つのは、孤立原子の被積分関数の対称性が高いためであって、一般の分子の場合にはFCC単独でこれほどの精度は出ないので、特にZの大きな原子を含む分子を扱う場合は、比較例２のｄｗ/ｄｖに数式２０や数式２３を使用する方法を併用する事が推奨される。

【0067】

【表1】

【実施例2】

【0068】

本実施例においては、主軸をZ方向に向けて配置されたSnTe分子のSn原子（下側）とTe原子（上側）に働く力を、両原子の核間距離をパラメーターとして計算した。DV-Xα法の標準の方式に従い、空間A_１と空間A_２において準乱数を用いて発生した標本点を数式６に従って空間Bに対応付けて７５２００個の標本点を得た。この分子軌道を用いて非特許文献３に記載された方法に従ってエネルギー勾配を計算した。こうして得られた結果を比較例として図３に示す。ここで核間距離３０１は実測の平衡核間距離（＝２．５２２８Å）を１．０とし、縦軸３０２は実施例１と同じ単位で表した両原子に働く力のｘ方向、y方向、ｚ方向の成分である。

【0069】

一方、空間A_１においてFCCの格子点を原点を含む３平面に対し面対称性を有する様に配列してSn原子用の標本点とし、同様に空間A_２においてFCCの格子点を原点を含む３平面に対し面対称性を有する様に配列してTe原子用の標本点とした。各標本点をR₀＝０．００２として数式２１に従って空間Bへ対応付けて得た７５２００個の標本点とその重みを用いて、DV-Xα法で分子軌道を求めた。この分子軌道を用いて非特許文献３に記載された方法に従ってエネルギー勾配を計算し、本実施例の結果として図３と同じ要領で図４に示す。

【0070】

図４においては、核間距離４０１が短いと下側に配置されたSn原子には下向きの力（４０７）が働き、上側に配置されたTe原子には絶対値が同じ上向きの力（４０８）が働き（斥力）、核間距離４０１が長いと下側に配置されたSn原子には上向きの力（４０７）が働き、上側に配置されたTe原子には絶対値が同じ下向きの力（４０８）が働く（引力）事が分かり、合理的な結果となっている。両原子に働く力がともに０となるのは核間距離＝１．００で、実測の平衡距離と良く一致した。一方図３においては、Sn原子に働く力のｚ成分３０７とTe原子に働く力のｚ成分３０８は定性的には図４と類似の傾向を示すものの、各核間距離３０１で両原子に働く力の絶対値は異なり、しかも本来存在しないはずのｘ方向の力（３０３，３０４）やｙ方向の力（３０５，３０６）がｚ方向の力（３０７や３０８）と同程度の大きさを持つと言う不合理な結果となった。本発明の数値積分法による精度向上の効果は明らかである。

【実施例3】

【0071】

本実施例においては、これまでに説明してきたエネルギー勾配を用いて、椅子型の６員環構造を有するSe_６分子の構造最適化を行った。使用した構造最適化の方法については非特許文献３に記載がある。まずSe_６分子に表２の欄１）に示すような平衡の構造と対称性は同じであるが、意図的にSe-Seの結合長を短く、∠SeSeSeを浅くした初期構造を与えた。次にDV-Xα法の標準の方式に従い、準乱数を用いて空間A_１〜空間A_６においてSe原子の標本点を生成した。各標本点を数式６に従って空間Bへ対応付けて得た１１２８００個の標本点とその重みを用いて、DV-Xα法で分子軌道を求めた。この分子軌道を用いて非特許文献３に記載された方法に従って各原子に働く力を計算した。各原子に働く力に合わせて各原子を少し移動させ、最初の工程と同様にして標本点と重みを得て分子軌道を求め各原子に働く力を再計算した。この様な工程を４回繰り返した所で、各原子に働く力に合わせて原子を動かしても全エネルギーが下がらなくなったので計算を終了した。最終の構造を表２の欄２）に示す。

【0072】

Se_６分子に対して比較例と同じ初期構造を与えた。空間A_１〜空間A_６においてFCCの格子点を生成して面対称性を有する標本点とし、これからR_０=０．００２として数式２１に従って空間Bへ対応付けて得た１１２８００個の標本点を用いて、DV-Xα法で分子軌道を求めた。この分子軌道を用いて非特許文献３に記載された方法に従って各原子に働く力を計算した。以下、比較例と同様にして構造最適化の工程を１６回繰り返した所で、各原子に働く力に合わせて原子を動かしても全エネルギーが下がらなくなったので計算を終了した。最終の構造を表２の欄３）に示す。あわせて実測によるSe_６分子の構造を表２の欄４）に示す。

【0073】

【表2】

【0074】

実施例の最終構造では、Se-Seの長さの誤差は-５％〜-２％以内で実測値とほぼ一致し、∠SeSeSeの誤差は-１．４°〜＋０．６°で実測値と一致した。一方比較例の最終構造では、Se-Seの長さの誤差は-４０％〜＋１７％に及び、∠SeSeSeの誤差は−７°〜＋３６°に及んだ。比較例のエネルギー勾配の計算精度が不十分なため、構造最適化の工程で適切な方向に構造を変化させる事が出来なかったためと思われる。本発明の方法による積分精度向上の効果は明らかである。

【実施例4】

【0075】

本発明の方法で標本点とその重みを生成し数値積分を行うプログラムをFortran言語で作成した。その概要を図５に示す。Mは被積分関数の特異点の数である。N_ｊは空間A_ｊに割り振る標本点の数であり、ｊに依存していてもよい。プログラムの５０１〜５１７で標本点とその重みを求める。計算の具体的な内容についてはこれまで説明してきた数式を引用して簡略に記述している。数式の番号の後の（ ）に別の番号がある場合、前の番号は空間Aの標本点ｐから空間Bへの標本点ｑへの変換を従来の数式６で行う場合の数式の番号を意味し、（ ）内の番号は標本点の変換を本発明の数式２１で行う場合の数式の番号を意味する。本実施例では被積分関数の分配比ｇ_ｉとして数式１４の関数を用いている。５１８〜５２４では得られた標本点とその重みを用いて数値積分を実行している。数値積分は全体を通して実行されるが、本発明による改良の本質的な部分は５０１〜５１７にある。しかし数値積分法として一体となったプログラムとする場合は、５２１を５１５の後に挿入しても良い。既存のプログラムで５１７の形で数値積分を行っている場合は、その前に５０1〜５１７をそのプログラムに挿入しても良い。図５から明らかなように、本実施例のプログラムは単純で計算の負荷は少ないと思われるが、大きなモデルの計算を行う等の理由で計算時間の短縮が求められる場合は、例えば５０６〜５１６のｄｏループをOpenMP等の手法により並列化して複数のコアを有する計算機にかければよい。

【0076】

実施例１から実施例３におけるエネルギー勾配の計算においては、５０１〜５１７行のプログラムで得た標本点とその重みをDV−Xα法のプログラムに読み込ませて所望の原子または分子の分子軌道を求めた。続いて非特許文献３に記載された方法でエネルギー勾配を求める既存のプログラムを全く改変なく使用して原子に働く力を求めた。既存の数値計算法を含むプログラムに本実施例のプログラムの結果を渡すだけで大幅な計算精度の向上が図れる点も本発明の特徴である。

【0077】

以上の実施例を通して、本発明の標本点とその重みを用いれば、値が大幅に変化する多変数の関数を、従来の方法よりはるかに精度良く積分出来る事を示した。一方従来の方法でも積分が可能だった関数でも、標本点の数が少なくとも同等の精度で積分できるので、計算時間を短縮する事も可能となる。

【産業上の利用可能性】

【0078】

クーロン力や万有引力等、空間で値が急激に変化する量の積分を伴う数値シミュレーションにおいて、これまで数値積分の誤差のため扱いが困難だったモデルが扱える、及び/又は必要な計算時間が大幅に短縮出来ると期待される。

【符号の説明】

【0079】

１０１ 面心立方格子の任意の一つの格子点（中心格子点）
１０２ 最近接格子点（１２個の内４個を示す。）
１０３ 中心格子点の空間要素（重み）
２０１ 標本点数
２０２ 誤差
２０３ 準乱数を用いて標本点を生成した場合の誤差
２０４ 面心立方格子を用いて標本点を生成した場合の誤差
２０１、３０１ 両原子の核間距離
２０２、３０２ 各原子に働く力
２０３、３０３ Sn原子に働く力のｘ成分
２０４、３０４ Te原子に働く力のｘ成分
２０５、３０５ Sn原子に働く力のｙ成分
２０６、３０６ Te原子に働く力のｙ成分
２０７、３０７ Sn原子に働く力のｚ成分
２０８、３０８ Te原子に働く力のｚ成分
５０１〜５１７ 本発明の標本点とその重みの生成方法を実施するためのプログラムのステートメント
５１８〜５２４ 本発明により与えられた標本点とその重みを用いて数値積分を行うプログラムのステートメント

【図1】

【図2】

【図3】

【図4】

【図5】