(57)【要約】【課題】一次相転移を回避しつつ負符号問題を解決することができるシミュレーション装置、コンピュータプログラム及びシミュレーション方法を提供する。【解決手段】複数のスピンの所定方向成分の和の平均を所定方向の磁化として演算する磁化演算部と、磁化の一次項及び二次以上の項を含む磁場関数を初期ハミルトニアンとして演算する初期ハミルトニアン演算部と、磁化とスピンの所定方向成分の和の平均との差を変数とするデルタ関数と磁場関数との乗算項を含む、初期ハミルトニアンに対する確率分布関数を演算する第１確率分布関数演算部と、スピン変数を更新するスピン変数更新部と、系が平衡状態になったか否かを判定する判定部と、系の平衡状態での所定方向の磁化を算出する第１磁化算出部と、所定方向の磁場を算出する磁場算出部と、磁場が定常状態であるか否かを判定する磁場判定部と、系に関する物理量を算出する物理量算出部とを備える。【選択図】図１
2018067200.tif 000002
【特許請求の範囲】
【請求項１】
二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートするシミュレーション装置であって、
前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する磁化演算部と、
該磁化演算部で演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する初期ハミルトニアン演算部と、
前記磁化演算部で演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する第１確率分布関数演算部と、
該第１確率分布関数演算部で演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新するスピン変数更新部と、
該スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する判定部と、
該判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する第１磁化算出部と、
該第１磁化算出部で算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する磁場算出部と、
該磁場算出部で算出した磁場が定常状態であるか否かを判定する磁場判定部と、
該磁場判定部で前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する物理量算出部と
を備えることを特徴とするシミュレーション装置。
【請求項２】
前記第１確率分布関数演算部で演算した確率分布関数に含まれる前記デルタ関数を積分表示にして、前記磁場関数の導関数を含む指数関数型演算子を用いて前記系のハミルトニアンに対する確率分布関数を演算する第２確率分布関数演算部を備え、
前記スピン変数更新部は、
前記第２確率分布関数演算部で演算して得られた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする請求項１に記載のシミュレーション装置。
【請求項３】
前記スピン変数更新部は、
前記磁場判定部で前記磁場が定常状態でないと判定した場合、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数を前記第１磁化算出部で算出した磁化に基づいて更新した前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする請求項２に記載のシミュレーション装置。
【請求項４】
予め前記所定方向の磁場の値を複数設定する設定部と、
該設定部で設定した磁場の値及び前記磁場関数の導関数の逆関数に基づいて前記所定方向の磁化を算出する第２磁化算出部と
を備え、
前記スピン変数更新部は、
前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数に前記設定部で設定した磁場の値を割り当てた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新し、
前記判定部は、
前記スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定し、
前記第１磁化算出部は、
前記判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出し、
前記物理量算出部は、
前記第１磁化算出部及び第２磁化算出部で算出した磁化が等しい場合、前記系に関する物理量を算出することを特徴とする請求項３に記載のシミュレーション装置。
【請求項５】
コンピュータに、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるコンピュータプログラムであって、
コンピュータに、
前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、
演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、
演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、
演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、
更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、
前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、
算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、
算出した磁場が定常状態であるか否かを判定する処理と、
前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理と
を実行させることを特徴とするコンピュータプログラム。
【請求項６】
二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるシミュレーション方法であって、
前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、
演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、
演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、
演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、
更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、
前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、
算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、
算出した磁場が定常状態であるか否かを判定する処理と、
前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理と
を含むことを特徴とするシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション装置、コンピュータプログラム及びシミュレーション方法に関する。
【背景技術】
【０００２】
離散的な多変数の１価関数（コスト関数）を最小化する問題を最適化問題という。パターン認識、自然言語処理、人工知能、機械学習をはじめとする大規模な計算を必要とする多くの重要な課題が最適化問題として定式化できる。また、近年、量子揺らぎなどの量子力学の性質を巧みに利用して最適化問題を解くアルゴリズムとしての量子アニーリングが注目されている。
【０００３】
量子アニーリングでは、コスト関数を二値変数の関数であるイジング（Ising）モデルとして表現し、その関数の最小値を求める問題として定式化する。量子アニーリングは、例えば、非特許文献１に記載されている。
【０００４】
量子アニーリングでは、式（１）で示すように、ｚ方向に向いたスピン（ベクトル表記法で示す太文字σで表すように多くのスピン）を反転させる量子力学的な揺らぎとして、ｘ方向の横磁場を利用する。この量子揺らぎの効果によって、より最適な解の探索を行う。
【０００５】
【数１】
2018067200.tif 000003
【０００６】
ここで、ハミルトニアンＨ０は、最適化問題のコスト関数を表現したものであり、ハミルトニアンＨ０の基底状態が最適解となるようにＨ０を選択する。σはスピン変数であり、スピンのｘ方向成分のシグマは、横磁場を表す初期ハミルトニアンであり、係数Γは、量子揺らぎの強さを制御するパラメータである。初期状態（時刻ｔ＝０）では、係数Γを非常に大きな値とし、時間の経過とともに係数Γを小さくし、最終的には０にする。最初は大きな量子揺らぎによって数多くの状態の重ね合わせを実現して状態探査をする。各時刻における瞬間的な基底状態を連続的にたどり、次第にΓが小さくなると、初期ハミルトニアンに比べてハミルトニアンＨ０の相対的な重みが大きくなり、最終的には、ハミルトニアンＨ０の基底状態に到達する。これは、最適化問題の解が得られ、所要の物理量が算出することができたということである。
【０００７】
図１５は量子系のエネルギーギャップの一例を示す模式図である。図１５において、横軸は時間を示し、縦軸はエネルギーを示す。量子系のエネルギー準位を時系列で追いかけた場合、基底状態と励起状態のエネルギー準位が接近する場合がある。Δは基底状態と第一励起状態とのエネルギーギャップを表す。Ｔは量子アニーリングに要する時間、すなわち最適解を求めるまでに要する計算時間を示す。また、Ｎはスピンの数を示す。図１５に示すようなエネルギーギャップΔでは、量子相転移（一次相転移）が起こり、計算時間Ｔは、指数関数的に増大して非常に大きな値となる。このため、量子アニーリングでも非常に計算時間が長いものも存在する。通常のコンピュータでも最適化問題を解くために非常に長い時間を要する問題があり、それらの問題に量子アニーリングを適用すると、同様に計算時間が長くなる。その背景には、この量子相移転の存在があるためである。
【０００８】
非特許文献１には、上述の一次相転移を回避する手法が記載されている。すなわち、式（２）に示すように、スピンのｘ方向成分の和の２乗の項（反強磁性ＸＸ相互作用とも称する）を追加する。式（２）において、ガンマγは係数である。
【０００９】
【数２】
2018067200.tif 000004
【００１０】
図１６は相図の一例を示す模式図である。図１６において、横軸は係数ガンマγであり、縦軸は係数Γの逆数である。符号ＱＰは量子常磁性相を示し、符号Ｆは強磁性相を示す。破線はＱＰとＦとの境界を示す。また、図１６中、横方向の実線は、一次相転移の線を示す。式（１）で示すハミルトニアンの場合には、符号Ａで示す破線のように、係数Γを非常に大きな値から小さい値にする場合、一次相転移の問題に遭遇する。一方、式（２）で示すハミルトニアンの場合には、符号Ｂで示す実線のように、一次相転移を回避することができる。
【先行技術文献】
【非特許文献】
【００１１】
【非特許文献１】Yuya Seki, Hidetoshi Nishimori, “Quantum annealing with antiferromagnetic fluctuations”, Phys. Rev. E 85, 051112(2012)
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかし、式（２）の反強磁性ＸＸ相互作用は、スピンのｘ方向成分の和の２乗の項を有しており、スピンのｘ方向成分による効果は、量子力学特有の効果であり、基本的には複素数で表現される。ＸＸ相互作用を利用した場合、その結果に複素数の２乗、すなわち負の値を生じる場合があり、確率的なサンプリングを行う際に必要なボルツマン重みが負となる、いわゆる負符号問題が生じ、通常のコンピュータ上でのシミュレーションが不可能とされていた。
【００１３】
本発明は斯かる事情に鑑みてなされたものであり、一次相転移を回避しつつ負符号問題を解決することができるシミュレーション装置、コンピュータプログラム及びシミュレーション方法を提供することを目的とする。
【課題を解決するための手段】
【００１４】
本発明の実施の形態に係るシミュレーション装置は、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートするシミュレーション装置であって、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する磁化演算部と、該磁化演算部で演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する初期ハミルトニアン演算部と、前記磁化演算部で演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する第１確率分布関数演算部と、該第１確率分布関数演算部で演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新するスピン変数更新部と、該スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する判定部と、該判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する第１磁化算出部と、該第１磁化算出部で算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する磁場算出部と、該磁場算出部で算出した磁場が定常状態であるか否かを判定する磁場判定部と、該磁場判定部で前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する物理量算出部とを備えることを特徴とする。
【００１５】
本発明の実施の形態に係るシミュレーション装置は、前記第１確率分布関数演算部で演算した確率分布関数に含まれる前記デルタ関数を積分表示にして、前記磁場関数の導関数を含む指数関数型演算子を用いて前記系のハミルトニアンに対する確率分布関数を演算する第２確率分布関数演算部を備え、前記スピン変数更新部は、前記第２確率分布関数演算部で演算して得られた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする。
【００１６】
本発明の実施の形態に係るシミュレーション装置は、前記スピン変数更新部は、前記磁場判定部で前記磁場が定常状態でないと判定した場合、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数を前記第１磁化算出部で算出した磁化に基づいて更新した前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする。
【００１７】
本発明の実施の形態に係るシミュレーション装置は、予め前記所定方向の磁場の値を複数設定する設定部と、該設定部で設定した磁場の値及び前記磁場関数の導関数の逆関数に基づいて前記所定方向の磁化を算出する第２磁化算出部とを備え、前記スピン変数更新部は、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数に前記設定部で設定した磁場の値を割り当てた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新し、前記判定部は、前記スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定し、前記第１磁化算出部は、前記判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出し、前記物理量算出部は、前記第１磁化算出部及び第２磁化算出部で算出した磁化が等しい場合、前記系に関する物理量を算出することを特徴とする。
【００１８】
本発明の実施の形態に係るコンピュータプログラムは、コンピュータに、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるコンピュータプログラムであって、コンピュータに、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを実行させることを特徴とする。
【００１９】
本発明の実施の形態に係るシミュレーション方法は、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるシミュレーション方法であって、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを含むことを特徴とする。
【発明の効果】
【００２０】
本発明によれば、一次相転移を回避しつつ負符号問題を解決することができる。
【図面の簡単な説明】
【００２１】
【図１】本実施の形態のシミュレーション装置の構成の一例を示す説明図である。
【図２】鈴木トロッター分解の一例を示す模式図である。
【図３】スピン変数の更新の様子の一例を示す模式図である。
【図４】横磁化ｍｘと横磁場ｍチルダｘとの関係の一例を示す説明図である。
【図５】本実施の形態のシミュレーション装置による適応的量子モンテカルロ法の処理手順の一例を示すフローチャートである。
【図６】本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第１例を示す説明図である。
【図７】本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第２例を示す説明図である。
【図８】本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第３例を示す説明図である。
【図９】本実施の形態のシミュレーション装置が行うデータ解析による量子モンテカルロ法の処理手順の一例を示すフローチャートである。
【図１０】本実施のデータ解析による量子モンテカルロ法の概念を示す説明図である。
【図１１】本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第１例を示す説明図である。
【図１２】本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第２例を示す説明図である。
【図１３】本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第３例を示す説明図である。
【図１４】本実施の形態のシミュレーション装置の構成の他の例を示す説明図である。
【図１５】量子系のエネルギーギャップの一例を示す模式図である。
【図１６】相図の一例を示す模式図である。
【発明を実施するための形態】
【００２２】
以下、本発明をその実施の形態を示す図面に基づいて説明する。図１は本実施の形態のシミュレーション装置１００の構成の一例を示す説明図である。本実施の形態のシミュレーション装置１００は、従来の量子アニーリングによるシミュレーションの範囲を大幅に広げることができ、イジング（Ising）モデルを用いて最適解を解くためのシミュレーションを実行することができる。
【００２３】
シミュレーション装置１００は、装置全体を制御する制御部１０、入力部１１、磁化演算部１２、初期ハミルトニアン演算部１３、第１確率分布関数演算部１４、第２確率分布関数演算部１５、スピン変数更新部１６、平衡状態判定部１７、出力部１８、第１磁化算出部１９、磁場算出部２０、磁場判定部２１、物理量算出部２２、第２磁化算出部２３、記憶部２４などを備える。
【００２４】
入力部１１は、シミュレーションを実行するための入力データ（例えば、トロッター数、スピンの数、横磁場の値など）を受け付ける。
【００２５】
出力部１８は、シミュレーションの結果である出力データ（例えば、エネルギーＥ、磁化ｍなど）を出力する。
【００２６】
磁化演算部１２は、複数のスピンの所定方向成分の和の平均を所定方向の磁化として演算する。ｚ方向に向いたスピンを反転させる量子力学的な揺らぎとして、ｘ方向の横磁場を利用する場合、所定方向は、横方向であるｘ方向とすることができ、所定方向の磁化は、横方向の磁化ｍｘとすることができる。以下では、所定方向を横方向として説明する。磁化演算部１２は、横磁化ｍｘを、スピンのｘ方向成分σｘのシグマの平均として演算する。
【００２７】
本実施の形態のシミュレーション装置１００は、式（３）において、右辺第２項で示すように、スピンのｘ方向成分σｘのシグマの平均を変数とする関数ｆにスピンの数Ｎを乗算したものを、初期ハミルトニアンとして定式化する。初期ハミルトニンアンは、ｚ方向に向いたスピンを反転させる量子力学的な揺らぎとして作用する。
【００２８】
【数３】
2018067200.tif 000005
【００２９】
具体的には、初期ハミルトニアン演算部１３は、横磁化ｍｘ一次項及び二次以上の項を含む磁場関数を初期ハミルトニアンとして演算する。式（４）に示すように、磁場関数をｆ（ｍｘ）で表すと、磁場関数ｆ（ｍｘ）は、例えば、ｆ（ｍｘ）＝Γ・ｍｘ＋（γ／２）・（ｍｘ）2 と表すことができる。ここで、係数Γは、量子揺らぎの強さを制御するパラメータであり、γは所定の係数である。初期ハミルトニアンは、磁場関数で表すことができる。磁場関数に横磁化ｍｘの２乗の項を含めることにより、一次相転移の問題を回避することが可能となる。なお、磁場関数は、式（４）に限定されるものではない。なお、上述の「二次以上の項」とは、二次項のみの場合、二次項に加えて三次項以上の高次項を含む場合、二次項は含まずに三次項以上の高次項を含む場合を意味する。
【００３０】
式（３）において、右辺の第１項である、ハミルトニアンＨ０は、最適化問題のコスト関数を表現したものであり、ハミルトニアンＨ０の基底状態が最適解となるようにＨ０を選択する。σはスピン変数であり、±１の値を取り得る。初期状態（時刻ｔ＝０）では、係数Γを非常に大きな値とし、時間の経過とともに係数Γを小さくし、最終的には０にする。最初は大きな量子揺らぎによって数多くの状態の重ね合わせを実現して状態探査をする。各時刻における瞬間的な基底状態を連続的にたどり、次第にΓが小さくなると、初期ハミルトニアンに比べてハミルトニアンＨ０の相対的な重みが大きくなり、最終的には、ハミルトニアンＨ０の基底状態に到達する。この状態で、最適化問題の解が得られ、所要の物理量を算出することができる。
【００３１】
第１確率分布関数演算部１４は、横磁化ｍｘとスピンの横方向成分σｘの和の平均との差を変数とするデルタ関数と、磁場関数との乗算項を含む指数関数型演算子を用いて初期ハミルトニアンに対する確率分布関数を演算する。具体的には、第１確率分布関数演算部１４によって演算される確率分布関数は、式（５）のように表すことができる。
【００３２】
【数４】
2018067200.tif 000006
【００３３】
式（５）は、鈴木トロッター分解したときの確率分布を示す。式（５）において、βは絶対温度の逆数に比例する係数であり、τはトロッター数を示す。式（５）に示すように、磁場関数ｆ（ｍｘ）に関する効果を書き直すことによって、ｆ（ｍｘ）を含む、任意の量子揺らぎ（式（３）の右辺第２項で表した初期ハミルトニアンの問題）を、横磁場を持つ単純なものに変更することができる。
【００３４】
式（５）に示すように、横磁化ｍｘとスピンのｘ方向成分σｘの和の平均との差を変数とするデルタ関数を導入することにより、横磁化ｍｘが、スピンのｘ方向成分σｘの和の平均と等しい場合には、デルタ関数が１となり、デルタ関数と磁場関数ｆ（ｍｘ）との乗算項は、結果として磁場関数ｆ（ｍｘ）に置き換えることができ、確率分布関数の指数関数型演算子からスピンのｘ方向成分（σｘ）の和の二次以上の高次項（ＸＸ相互作用を含む）を除くことができるので、負符号問題がなくなる。
【００３５】
式（６）は、鈴木トロッター分解公式を示す。式（６）において、Ａ、Ｂは演算子であり、Ｌはトロッター数である。なお、本明細書では、トロッター数をτでも表している。
【００３６】
【数５】
2018067200.tif 000007
【００３７】
一般的に、量子系のハミルトニアンは、局所的な構成要素間の相互作用を表す局所ハミルトニアンの和で定義することができる。局所ハミルトニアン同士はお互いに非可換であるため、量子系のハミルトニアンの行列表現のサイズが大きくなり、計算コストが高くなる。そこで、式（６）で示す鈴木トロッター分解公式を用いて、指数演算子を行列表現のサイズが小さな局所ハミルトニアンの指数演算子の積に分解することができる。
【００３８】
第２確率分布関数演算部１５は、第１確率分布関数演算部１４で演算した確率分布関数に含まれるデルタ関数を積分表示にして、磁場関数の導関数を含む指数関数型演算子を用いて系のハミルトニアンに対する確率分布関数を演算する。具体的には、第２確率分布関数演算部１５によって演算される確率分布関数は、式（７）のように表すことができる。
【００３９】
【数６】
2018067200.tif 000008
【００４０】
式（７）において、τはトロッター数であり、ｍチルダｘは、式（８）のように、磁場関数ｆ（ｍｘ）のｍｘについての導関数ｆ′（ｍｘ）で表すことができる。式（７）は、式（５）のデルタ関数を積分表示（フーリエ積分表示）することで、ｍチルダｘが導入されている。式（７）では、スピンのｘ方向成分σｘを持つハミルトニアンの問題を横磁場の問題として扱うことができ、鈴木トロッター分解を実行することによって、スピンのｘ方向成分σｘに関する項を式（９）のＢを係数として持つトロッター間相互作用に置き換えることができ、数値計算を実行することが可能となる。これにより、通常のコンピュータを用いてシミュレーションを実行することができる。なお、係数Ｂは、式（１０）で表すことができる。また、本明細書では、式（９）も第２確率分布関数演算部１５によって演算される確率分布関数と称する。式（９）において、Ｚ（ｍｘ）は規格化のための係数である。
【００４１】
図２は鈴木トロッター分解の一例を示す模式図である。図２において、横軸は１次元格子上でのスピン変数を表し、いわゆる実空間方向を示す。縦軸は鈴木トロッター分解によって導入された方向（トロッター方向）であり、２次元の格子点上に状態変数が配置される。このように、鈴木トロッター分解によって量子モデルは、次元が一つ増えた状態空間を持つ古典モデルに変換されたと考えることができる。
【００４２】
式（７）の指数関数部分に注目して、ｍｘについて極大値を求める。ｍｘについての極大を取るところだけに注目する鞍点法を利用することで、式（９）を導出することができる。その鞍点の条件（極大点となる）を表す式が、式（８）となり、式（９）でｍチルダｘが消えている。式（９）は、横磁化ｍｘが決まるとスピン変数σがどのような振る舞いをするかを示す。
【００４３】
スピン変数更新部１６は、第１確率分布関数演算部１４によって演算して得られた確率分布に基づいて、複数のスピンそれぞれのスピン変数を更新する。より具体的は、スピン変数更新部１６は、第２確率分布関数演算部１５で演算して得られた系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。スピン変数の更新は、例えば、スピン変数を選び、熱浴法又はメトロポリス法などの詳細釣り合いを満たすもので更新することができる。
【００４４】
図３はスピン変数の更新の様子の一例を示す模式図である。左側の図は、現在のスピン変数の状態を示し、右図は任意のスピン変数として、格子ｉのスピン変数を選んだとする。なお、便宜上、格子ｉのスピン変数の周囲に４つのスピン変数があるとし、格子ｉのスピン変数しか変化できないとする（１スピンフリップ）。
【００４５】
まず、左図の状態でのスピン変数を式（９）に代入して、現在の確率分布Ｐｐを算出する。次に、右図のように、選択した格子ｉのスピン変数を変更した状態でのスピン変数を式（９）に代入して、次の新しい確率分布Ｐｎを算出する。フリップ確率ＰｆをＰｆ＝Ｐｎ／Ｐｐにより算出する。一様乱数ｒを生成し、フリップ確率Ｐｆが乱数ｒより大きい場合、格子ｉのスピン変数を右図のとおり（スピンをフリップさせる）とし、フリップ確率Ｐｆが乱数ｒより大きくない場合、格子ｉのスピン変数を左図のとおり（スピンをフリップさせない）とする。上述の処理をすべての格子（すなわち、スピンの数Ｎ）について繰り返す。
【００４６】
平衡状態判定部１７は、判定部としての機能を有し、更新したスピン変数に基づいて系が平衡状態になったか否かを判定する。系が平衡状態になったか否かの判定は、例えば、系のエネルギーＥ、磁化ｍ（磁化の２乗、４乗などのｍ−メント）を算出（測定）して、変動がなければ系が平衡状態になったと判定することができる。系のエネルギーＥは、式（１）で算出する。磁化ｍは、全スピンについて和を計算し、その和をスピンの数で除算して求めることができる。
【００４７】
第１磁化算出部１９は、平衡状態判定部１７で系が平衡状態になったと判定した場合、当該平衡状態での横磁化ｍｘを算出する。平衡状態での横磁化ｍｘは、式（１１）により算出した量の時間平均値で求めることができる。式（１１）において、ｉはスピンの場所（格子点）を示し、ｔはトロッター数を示し、τはトロッター総数を示し、Ｎはスピン総数を示す。
【００４８】
【数７】
2018067200.tif 000009
【００４９】
磁場算出部２０は、第１磁化算出部１９で算出した横磁化ｍｘ及び磁場関数ｆ（ｍｘ）に基づいて、複数のスピンに対する横磁場を算出する。横磁場（ｍチルダｘ）は、磁場関数ｆ（ｍｘ）を横磁化ｍｘで微分した式（８）に横磁化ｍｘを代入して算出することができる。
【００５０】
磁場判定部２１は、磁場算出部２０で算出した横磁場（ｍチルダｘ）が定常状態であるか否かを判定する。横磁場が定常状態であるか否かは、平衡状態での横磁化ｍｘの値に応じて横磁場の値を適応的に変化させ、横磁場が変化しない場合、定常状態であると判定することができる。
【００５１】
物理量算出部２２は、磁場判定部２１で横磁場が定常状態であると判定した場合、系に関する物理量を算出する。系に関する物理量は、例えば、エネルギーＥ、磁化ｍなどである。シミュレーションを続けながら、スピン変数に基づいて、エネルギーＥ、磁化ｍなどの物理量を繰り返し計算し、ある程度の時間が経過したときに物理量の時間平均を算出し、物理量の結果とする。時間平均によって任意の精度で計算することができ、時間を長くすれば、それに応じて精度を高めることができる。
【００５２】
記憶部２４は、シミュレーションの必要なデータ、入力データ、シミュレーション中に得られた処理結果、出力データなどを記憶する。
【００５３】
上述の構成により、一次相転移を回避しつつ負符号問題を解決して、シミュレーションを実施することができる。
【００５４】
また、スピン変数更新部１６は、磁場判定部２１で横磁場が定常状態でないと判定した場合、系のハミルトニアンに対する確率分布関数に含まれる磁場関数ｆ（ｍｘ）の導関数ｆ′（ｍｘ）を第１磁化算出部１９で算出した横磁化に基づいて更新し、更新した系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。
【００５５】
具体的には、式（９）の中にある係数Ｂが、式（１０）で示すように、磁場関数ｆ（ｍｘ）の導関数ｆ′（ｍｘ）を含んでいる。すなわち、ｍｘによってｆ′（ｍｘ）を更新すると、係数Ｂが変わり、結果として式（９）で演算される系の確率分布が変わる。
【００５６】
図４は横磁化ｍｘと横磁場ｍチルダｘとの関係の一例を示す説明図である。ｍチルダｘは、磁場関数ｆ（ｍｘ）の導関数ｆ′（ｍｘ）である。例えば、磁場関数ｆ（ｍｘ）が式（４）で表される場合、導関数ｆ′（ｍｘ）は、横磁化ｍｘの一次関数となる。すなわち、横磁化ｍｘを変えることによって横磁場（ｍチルダｘ）を変えることができる。なお、従来の量子アニーリングでは、横磁場は一定であった。
【００５７】
上述のように、横磁場が定常状態でない場合、解ではない別の横磁場又は横磁化が得られたことになる。そこで、系のハミルトニアンに対する確率分布を横磁化に応じて計算し、再度スピン変数を更新する処理を行うことにより、解に到達することができる。
【００５８】
次に、本実施の形態のシミュレーション装置１００の動作について説明する。図５は本実施の形態のシミュレーション装置１００による適応的量子モンテカルロ法の処理手順の一例を示すフローチャートである。適応的量子モンテカルロ法での「適応的」とは、従来の量子アニーリング法を実現するための確率的手法である量子モンテカルロ法と区別するための用語であり、量子相移転（一次相転移）を回避しつつ負符号問題も回避することができ、通常のコンピュータにおいて実行可能なシミュレーション法であることを表す。なお、以下では、便宜上、処理の主体を制御部１０として説明する。
【００５９】
制御部１０は、トロッター数、スピンの数を設定する（Ｓ１１）。トロッター数は、シミュレータ（コンピュータ）の性能にも依存するが、例えば、１２８とすることができる。スピンの数Ｎは、任意のサイズとすることができる。トロッター数もスピンの数も大きくすれば、計算精度を高めることができる。
【００６０】
制御部１０は、スピン変数のシグマに初期値を設定する（Ｓ１２）。スピン変数のシグマに初期値を設定すると、横磁化ｍｘを計算することができるので、横磁場の値（ｍチルダｘ、すなわちｆ′（ｍｘ））を初期値に設定することができる。
【００６１】
制御部１０は、スピン変数を乱数によって初期化し（Ｓ１３）、スピン変数を選び、熱浴法又はメトロポリス法（詳細釣り合いを満たすもの）で更新する（Ｓ１４）。これにより、系のハミルトニアンが目的ハミルトニアンＨ０の基底状態に向かって収束するようにする。
【００６２】
制御部１０は、系が平衡状態になったか否かを判定する（Ｓ１５）。平衡状態になったか否かの判定は、エネルギーＥの値、磁化ｍなどが変動するか否かによって行うことができる。
【００６３】
系が平衡状態でない場合（Ｓ１５でＮＯ）、制御部１０は、ステップＳ１４以降の処理を続ける。系が平衡状態である場合（Ｓ１５でＹＥＳ）、制御部１０は、平衡状態での横磁化ｍｘの値に応じて横磁場（ｍチルダｘ、すなわちｆ′（ｍｘ））の値を変更する（Ｓ１６）。
【００６４】
制御部１０は、横磁場が変化するか、すなわち横磁場が定常状態であるか否かを判定する（Ｓ１７）。横磁場が変化する場合（Ｓ１７でＹＥＳ）、制御部１０は、解ではない別の横磁場又は横磁化が得られたとして、ステップＳ１４以降の処理を続ける。
【００６５】
横磁場が変化しない場合（Ｓ１７でＮＯ）、制御部１０は、解が得られたとして物理量の測定（計算）を開始し（Ｓ１８）、所要の時間が経過したときに、測定量の時間平均を求め、物理量の結果とし、処理を終了する。
【００６６】
図６は本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第１例を示す説明図である。図６において、横軸は横磁場Γ（Gamma）を示し、縦軸はエネルギーＥを示す。図中Ｎはシミュレーション時のスピンの数を示し、exactγ＝１で示す実線は、手計算により求めた真の解（正解）を示す。図６から分かるように、シミュレーション結果は、真の解とほぼ一致し、特にスピン数Ｎを大きくすることにより、より真の解に近づくことが分かる。
【００６７】
図７は本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第２例を示す説明図であり、図８は本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第３例を示す説明図である。図７において、縦軸は磁化ｍを示す。また、図８において、縦軸は横磁化ｍｘを示す。図７及び図８においても、シミュレーション結果は、スピン数Ｎを大きくすることにより、より真の解に近づくことが分かる。
【００６８】
次に、本実施の形態のシミュレーション装置１００が実行するテータ解析による量子モンテカルロ法について説明する。テータ解析による量子モンテカルロ法は、通常の量子モンテカルロ法である。
【００６９】
制御部１０は、設定部としての機能を有し、予め横磁場の値を複数設定する。すなわち、横磁場の値を複数用意しておく。
【００７０】
第２磁化算出部２３は、横磁場の値及び磁場関数の導関数の逆関数に基づいて横磁化を算出する。例えば、横磁場ｍチルダｘ＝ｆ′（ｍｘ）とすると、横磁化ｍｘは、ｍｘ＝ｆ′（ｍチルダｘ）の逆関数で算出することができる。複数の横磁場に対して横磁化を算出することにより、横磁場と横磁化との関係をプロットすることができる。
【００７１】
スピン変数更新部１６は、系のハミルトニアンに対する確率分布関数に含まれる磁場関数の導関数に制御部１０が設定した横磁場の値を割り当てた系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。
【００７２】
平衡状態判定部１７は、スピン変数更新部１６で更新したスピン変数に基づいて系が平衡状態になったか否かを判定する。
【００７３】
第１磁化算出部１９は、系が平衡状態になった場合、当該平衡状態での横磁化を算出する。すなわち、予め容易した横磁場を用いて量子モンテカルロ法を実行して、横磁化を算出する。これにより、横磁場と横磁化との関係をプロットすることができる。
【００７４】
物理量算出部２２は、第１磁化算出部１９及び第２磁化算出部２３で算出した横磁化が等しい場合、系に関する物理量を算出する。すなわち、磁場関数の導関数の逆関数に基づいて算出した横磁化と、量子モンテカルロ法を実行して計算された横磁化とが等しくなる場合、当該横磁化及び対応する横磁場を求めることにより、データ解析的なアプローチによっても解を求めることができる。
【００７５】
図９は本実施の形態のシミュレーション装置１００が行うデータ解析による量子モンテカルロ法の処理手順の一例を示すフローチャートである。制御部１０は、横磁場（ｍチルダｘ、すなわち、ｆ′（ｍｘ））の値を複数設定し（Ｓ３１）、設定したそれぞれの横磁場の値を用いて量子モンテカルロ法を実行する（Ｓ３２）。
【００７６】
制御部１０は、量子モンテカルロ法を実行することによって得られた横磁化の値と、当該横磁化が得られたときの横磁場とを対応付けて横磁場の値と横磁化の値との関係をプロットする（Ｓ３３）。なお、ここで、プロットするとは、実際にチャートに描く必要はなく、対応付けが分かるような形態であればよい。
【００７７】
制御部１０は、磁場関数ｆ（ｍｘ）の導関数ｆ′（ｍｘ）の逆関数に基づいて横磁化を算出する（Ｓ３４）。具体的には、設定した横磁場の値それぞれを逆関数に代入して横磁化を算出することができる。
【００７８】
制御部１０は、逆関数によって算出した横磁化が、プロット上の横磁化と一致するときの横磁化の値及び横磁場の値を特定する（Ｓ３５）。これにより、真の解が得られ、適応的量子モンテカルロ法と同様に物理量の結果を得ることができる。制御部１０は、処理を終了する。
【００７９】
図１０は本実施のデータ解析による量子モンテカルロ法の概念を示す説明図である。図１０において、横軸は横磁場（ｍチルダｘ）の値を示し、縦軸は横磁化ｍｘの値を示す。図１０中、符号Ｐ１で示す曲線は、量子モンテカルロ法によって得られた結果をプロットしたものである。また、横磁場ｍチルダｘ＝ｆ′（ｍｘ）とすると、横磁化ｍｘは、ｍｘ＝ｆ′（ｍチルダｘ）の逆関数で算出することができる。符号Ｐ２で示す直線は、ｍｘ＝ｆ′（ｍチルダｘ）の逆関数を満たす。そして、曲線Ｐ１と直線Ｐ２とが交差する点における横磁化ｍｘの値と横磁場ｍチルダｘの値が解となる。
【００８０】
図１１は本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第１例を示す説明図である。図１１において、横軸は横磁場Γ（Gamma）を示し、縦軸はエネルギーＥを示す。図中Ｎはシミュレーション時のスピンの数を示し、各Ｎに対応するチャートは、曲線Ｐ１と直線Ｐ２とが交差する点の集合である。exactγ＝１で示す実線は、手計算により求めた真の解（正解）を示す。図１１から分かるように、シミュレーション結果は、真の解とほぼ一致し、特にスピン数Ｎを大きくすることにより、より真の解に近づくことが分かる。
【００８１】
図１２は本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第２例を示す説明図であり、図１３は本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第３例を示す説明図である。図１２において、縦軸は磁化ｍを示す。また、図１３において、縦軸は横磁化ｍｘを示す。図１２及び図１３においても、シミュレーション結果は、スピン数Ｎを大きくすることにより、より真の解に近づくことが分かる。
【００８２】
本実施の形態のシミュレーション装置１００は、ＣＰＵ（プロセッサ）、ＲＡＭ（メモリ）などを備えたコンピュータを用いて実現することもできる。すなわち、図５及び図９に示すような、各処理の手順を定めたコンピュータプログラムをコンピュータに備えられたＲＡＭ（メモリ）にロードし、コンピュータプログラムをＣＰＵ（プロセッサ）で実行することにより、シミュレーション装置１００を実現することができる。
【００８３】
図１４は本実施の形態のシミュレーション装置の構成の他の例を示す説明図である。図１４において、符号３００は、通常のコンピュータである。コンピュータ３００は、制御部３０、入力部４０、出力部５０、外部Ｉ／Ｆ（インタフェース）部６０などを備える。制御部３０は、ＣＰＵ３１、ＲＯＭ３２、ＲＡＭ３３、Ｉ／Ｆ（インタフェース）３４などを備える。
【００８４】
入力部４０は、シミュレーションのための入力データを取得する。出力部５０は、シミュレーション結果である出力データを出力する。Ｉ／Ｆ３４は、制御部３０と、入力部４０、出力部５０及び外部Ｉ／Ｆ部６０それぞれとの間のインタフェース機能を有する。
【００８５】
外部Ｉ／Ｆ部６０は、コンピュータプログラムを記録した記録媒体Ｍ（例えば、ＤＶＤなどのメディア）からコンピュータプログラムを読み取ることが可能である。
【００８６】
なお、図示していないが、記録媒体Ｍに記録されたコンピュータプログラムは、持ち運びが自由なメディアに記録されたものに限定されるものではなく、インターネット又は他の通信回線を通じて伝送されるコンピュータプログラムも含めることができる。また、コンピュータには、複数のプロセッサを搭載した１台のコンピュータ、あるいは、通信ネットワークを介して接続された複数台のコンピュータで構成されるコンピュータシステムも含まれる。
【００８７】
上述のように、本実施の形態によれば、シミュレーションの計算時間が最適解を得ることができないほど長くなる量子相転移の問題を回避することができる。また、本実施の形態によれば、確率密度が負になってしまうという負符号問題も回避することができるので、限定的なモデルに限られることなく、幅広い範囲のモデルに対して、通常のコンピュータ上でシミュレーションを実行することが可能となり、量子モンテカルロ法の適用範囲を広げることができ、物質探索、あるいは設計のシミュレーションの範囲を拡大させることができる。また、本実施の形態のシミュレーション方法は、量子コンピュータの製作現場、人工知能、機械学習などの大規模な計算を必要とする技術分野で利用することができる。
【００８８】
本実施の形態に係るシミュレーション装置は、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートするシミュレーション装置であって、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する磁化演算部と、該磁化演算部で演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する初期ハミルトニアン演算部と、前記磁化演算部で演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する第１確率分布関数演算部と、該第１確率分布関数演算部で演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新するスピン変数更新部と、該スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する判定部と、該判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する第１磁化算出部と、該第１磁化算出部で算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する磁場算出部と、該磁場算出部で算出した磁場が定常状態であるか否かを判定する磁場判定部と、該磁場判定部で前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する物理量算出部を備えることを特徴とする。
【００８９】
本実施の形態に係るコンピュータプログラムは、コンピュータに、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるコンピュータプログラムであって、コンピュータに、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを実行させることを特徴とする。
【００９０】
本実施の形態に係るコンピュータでの読み取り可能な記録媒体は、コンピュータに、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるコンピュータプログラムであって、コンピュータに、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを実行させるためのコンピュータプログラムを記録したことを特徴とする。
【００９１】
本実施の形態に係るシミュレーション方法は、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるシミュレーション方法であって、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを含むことを特徴とする。
【００９２】
磁化演算部は、複数のスピンの所定方向成分の和の平均を当該所定方向の磁化として演算する。ｚ方向に向いたスピンを反転させる量子力学的な揺らぎとして、ｘ方向の横磁場を利用する場合、所定方向は、横方向であるｘ方向とすることができ、所定方向の磁化は、横方向の磁化ｍｘとすることができる。すなわち、横磁化ｍｘを、スピンのｘ方向成分のシグマの平均として演算する。
【００９３】
初期ハミルトニアン演算部は、演算した磁化の一次項及び二次以上の項を含む磁場関数を初期ハミルトニアンとして演算する。磁場関数をｆ（ｍｘ）で表すと、磁場関数ｆ（ｍｘ）は、例えば、ｆ（ｍｘ）＝Γ・ｍｘ＋（γ／２）・（ｍｘ）2 と表すことができる。初期ハミルトニアンは、磁場関数で表すことができる。磁場関数に横磁化ｍｘの２乗の項を含めることにより、一次相転移の問題を回避することが可能となる。
【００９４】
第１確率分布関数演算部は、演算した磁化とスピンの所定方向成分の和の平均との差を変数とするデルタ関数と、磁場関数との乗算項を含む指数関数型演算子を用いて初期ハミルトニアンに対する確率分布関数を演算する。磁化とスピンの所定方向成分の和の平均との差を変数とするデルタ関数を導入することにより、横磁化ｍｘが、スピンの所定方向成分の和の平均と等しい場合には、デルタ関数が１となり、デルタ関数と磁場関数との乗算項は、結果として磁場関数に置き換えられるので、確率分布関数の指数関数型演算子からスピンの所定方向成分（σｘ）の和の２次以上の高次項を除くことができるので、負符号問題がなくなる。
【００９５】
スピン変数更新部は、演算して得られた確率分布に基づいて、複数のスピンそれぞれのスピン変数を更新する。例えば、スピン変数を選び、熱浴法又はメトロポリス法などの詳細釣り合いを満たすもので更新することができる。
【００９６】
判定部は、更新したスピン変数に基づいて系が平衡状態になったか否かを判定する。系が平衡状態になったか否かの判定は、例えば、系のエネルギー、磁化（磁化の２乗、４乗などのｍ−メント）を算出（測定）して、変動がなければ系が平衡状態になったと判定することができる。
【００９７】
第１磁化算出部は、系が平衡状態になった場合、当該平衡状態での所定方向の磁化（横磁化ｍｘ）を算出する。
【００９８】
磁場算出部は、算出した磁化（横磁化）及び磁場関数に基づいて、複数のスピンに対する所定方向の磁場（横磁場）を算出する。横磁場は、磁場関数を横磁化で微分した式に横磁化を代入して算出することができる。
【００９９】
磁場判定部は、算出した磁場（横磁場）が定常状態であるか否かを判定する。横磁場が定常状態であるか否かは、平衡状態での横磁化の値に応じて横磁場の値を適応的に変化させ、横磁場が変化しない場合、定常状態であると判定することができる。
【０１００】
物理量算出部は、横磁場が定常状態であると判定した場合、系に関する物理量を算出する。系に関する物理量は、例えば、エネルギー、磁化などである。シミュレーションを続けながら、スピン変数に基づいて、エネルギー、磁化などの物理量を繰り返し計算し、ある程度の時間が経過したときに物理量の時間平均を算出し、物理量の結果とする。
【０１０１】
上述の構成により、一次相転移を回避しつつ負符号問題を解決して、シミュレーションを実施することができる。
【０１０２】
本実施の形態に係るシミュレーション装置は、前記第１確率分布関数演算部で演算した確率分布関数に含まれる前記デルタ関数を積分表示にして、前記磁場関数の導関数を含む指数関数型演算子を用いて前記系のハミルトニアンに対する確率分布関数を演算する第２確率分布関数演算部を備え、前記スピン変数更新部は、前記第２確率分布関数演算部で演算して得られた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする。
【０１０３】
第２確率分布関数演算部は、第１確率分布関数演算部で演算した確率分布関数に含まれるデルタ関数を積分表示にして、磁場関数の導関数を含む指数関数型演算子を用いて系のハミルトニアンに対する確率分布関数を演算する。デルタ関数を積分表示にすることにより、磁場関数の導関数を導入することができ、鈴木トロッター分解を実行することで、スピンのｘ方向成分（σｘ）に関する項をトロッター間相互作用に置き換えることができ、数値計算を実行することが可能となる。
【０１０４】
スピン変数更新部は、第２確率分布関数演算部で演算して得られた系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。これにより、通常のコンピュータを用いてシミュレーションを実行することができる。
【０１０５】
本実施の形態に係るシミュレーション装置は、前記スピン変数更新部は、前記磁場判定部で前記磁場が定常状態でないと判定した場合、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数を前記第１磁化算出部で算出した磁化に基づいて更新した前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする。
【０１０６】
スピン変数更新部は、横磁場が定常状態でないと判定した場合、系のハミルトニアンに対する確率分布関数に含まれる磁場関数の導関数を第１磁化算出部で算出した横磁化に基づいて更新した系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。
【０１０７】
横磁場が定常状態でない場合、解ではない別の横磁場又は横磁化が得られたことになる。そこで、系のハミルトニアンに対する確率分布を横磁化に応じて計算し、再度スピン変数を更新する処理を行うことにより、解に到達することができる。
【０１０８】
本実施の形態に係るシミュレーション装置は、予め前記所定方向の磁場の値を複数設定する設定部と、該設定部で設定した磁場の値及び前記磁場関数の導関数の逆関数に基づいて前記所定方向の磁化を算出する第２磁化算出部とを備え、前記スピン変数更新部は、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数に前記設定部で設定した磁場の値を割り当てた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新し、前記判定部は、前記スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定し、前記第１磁化算出部は、前記判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出し、前記物理量算出部は、前記第１磁化算出部及び第２磁化算出部で算出した磁化が等しい場合、前記系に関する物理量を算出することを特徴とする。
【０１０９】
設定部は、予め所定方向の磁場（横磁場）の値を複数設定する。すなわち、横磁場の値を複数用意しておく。
【０１１０】
第２磁化算出部は、横磁場の値及び磁場関数の導関数の逆関数に基づいて所定方向の磁化（横磁化）を算出する。複数の横磁場に対して横磁化を算出することにより、横磁場と横磁化との関係をプロットすることができる。
【０１１１】
スピン変数更新部は、系のハミルトニアンに対する確率分布関数に含まれる磁場関数の導関数に設定部で設定した磁場の値を割り当てた系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。判定部は、スピン変数更新部で更新したスピン変数に基づいて系が平衡状態になったか否かを判定する。第１磁化算出部は、系が平衡状態になった場合、当該平衡状態での横磁化を算出する。すなわち、予め容易した横磁場を用いて量子モンテカルロ法を実行して、横磁化を算出する。これにより、横磁場と横磁化との関係をプロットすることができる。
【０１１２】
物理量算出部は、第１磁化算出部及び第２磁化算出部で算出した磁化が等しい場合、系に関する物理量を算出する。すなわち、磁場関数の導関数の逆関数に基づいて算出した横磁化と、量子モンテカルロ法を実行して計算された横磁化とが等しくなる場合、当該横磁化及び対応する横磁場を求めることにより、データ解析的なアプローチによっても解を求めることができる。
【符号の説明】
【０１１３】
１０、３０  制御部
１１、４０  入力部
１２  磁化演算部
１３  初期ハミルトニアン演算部
１４  第１確率分布関数演算部
１５  第２確率分布関数演算部
１６  スピン変数更新部
１７  平衡状態判定部
１８、５０  出力部
１９  第１磁化算出部
２０  磁場算出部
２１  磁場判定部
２２  物理量算出部
２３  第２磁化算出部
２４  記憶部
３１  ＣＰＵ
３２  ＲＯＭ
３３  ＲＡＭ
３４  Ｉ／Ｆ
６０  外部Ｉ／Ｆ部
【図１】
2018067200.tif 000010
【図２】
2018067200.tif 000011
【図３】
2018067200.tif 000012
【図４】
2018067200.tif 000013
【図５】
2018067200.tif 000014
【図６】
2018067200.tif 000015
【図７】
2018067200.tif 000016
【図８】
2018067200.tif 000017
【図９】
2018067200.tif 000018
【図１０】
2018067200.tif 000019
【図１１】
2018067200.tif 000020
【図１２】
2018067200.tif 000021
【図１３】
2018067200.tif 000022
【図１４】
2018067200.tif 000023
【図１５】
2018067200.tif 000024
【図１６】
2018067200.tif 000025