原子の「谷」はどこにある？――H100×2でNbSe2の388構造を下見し、DFT 5点で確かめた

結論から

今回調べたのは、「H100なら計算が速い」という話だけではありません。調和フォノンとfrozen-modeの数値は、論文とともに公開された振動特性向け機械学習ポテンシャル「V1」上の結果です。Quantum ESPRESSOによるDFT――電子の状態から原子に働く力とエネルギーを計算する、より精密な方法――の結果は、後半で明示して分けます。

二セレン化ニオブ（NbSe2）の原子が、電荷密度波（CDW）と呼ばれる周期構造へ向かって本当に動きたがっているのかを、計算の大きさと条件を変えながら確かめました。

今回使った「V1」とは

ここでいうV1は、H100の機能名でも、ソフトウェアのバージョン番号でもありません。公開データでV1_monolayer_vibrationalと識別されている、単層NbSe2の原子振動を調べるための学習済みAIモデルです。原子シミュレーション向けニューラルネットワークAllegroを使い、第一原理計算のデータから「原子の並び方」と「そのときのエネルギー・各原子に働く力」の関係を学んでいます。

文章や画像を作る生成AIではありません。原子配置を入力すると、エネルギーと力を返す、材料計算専用の予測器です。DFTが一つ一つの配置について電子状態から計算するのに対し、V1は学習済みの関係を使って高速に予測します。そのため多数の候補を調べられますが、DFTそのものではなく、学習データとモデルが表現できる範囲に基づく近似です。未知の全構造での正しさや、最終的なCDW構造をV1だけで保証することはできません。

今回の計算で、各ソフトは何をしたのか

• AllegroとV1：Allegroは原子の周囲の情報からエネルギーと力を予測するモデルの仕組みで、V1はその仕組みを使ってNbSe2の振動特性向けに学習した完成済みモデルです。今回、H100が主に高速化したのは、このV1による多数のエネルギー・力計算です。
• ASEとPhonopy：ASEは原子構造を作ってV1へ渡し、結果を受け取る進行役です。Phonopyは原子をわずかに動かした構造を用意し、V1が返した力からフォノン、つまり原子集団の振動を組み立てる解析役です。388構造の探索とフォノン計算は、このPython/ASE経路が中心です。
• LAMMPS（ランプス）：Allegroモデルを組み込んでMDなどを実行できる原子シミュレーションソフトです。今回はMDの主計算には使わず、代表構造で同じV1を別経路から動かし、Python/ASE経路とエネルギー・力が一致するかを確かめる実装検査に使いました。
• Quantum ESPRESSO（QE）：電子状態からエネルギーと力を計算する第一原理DFTソフトです。H100ではなく別ホストのCPU 64 MPI ranksで動かし、後半の独立確認では、V1で選んだ27原子の5点を精密に計算し直しました。

流れを一行で表すと、ASE・Phonopyが検査条件を作る → Allegro/V1がH100で多数の候補を評価する → Quantum ESPRESSOが重要な少数点をDFTで確かめる、となります。QEをH100で高速化した計算でも、LAMMPSの結果をQEへ変換した計算でもありません。

H100がなければ原理的に不可能という意味ではありません。今回H100×2がもたらした価値は、異なる構造や条件を2枚へ分担し、多数の候補と不合格条件まで現実的な時間で検査できたことです。

この計算を一言でいうと

材料の中で原子が取りうる並び方は、山や谷が無数にある広い地形に似ています。今回の目的は、最終的な谷底をいきなり決めることではなく、原子をどちらへ少しずらすとエネルギーが下がり、選んだ一本の道筋のどこで坂の向きが変わるかを見つけることです。計算箱を固定し、CDWの一つの振付と開始位置だけをたどる静的な検査であり、原子の自然な時間変化や、あらゆる構造候補を調べたものではありません。

• 知りたいこと：選んだCDWの振付に沿って、少しずらすとエネルギーが下がる方向は本当にあるか
• 難しいこと：計算箱の大きさ、動かす幅、波の開始位置を変えると、候補数が掛け算で増える
• 今回の方法：H100上のV1では3×3～12×12、振幅、位相、対照を含む388構造を検査。これと同じ仕事の速度比較ではなく、選んだ27原子の一次元経路5点をDFTで計算し直す
• 得られたこと：正負両方向で、経路の傾きがゼロになる地点を含む範囲を挟めた。ただし、全方向で最も安定な構造はまだ決めていない

これが実用上重要なのは、次の高精度な構造緩和や有限温度計算を、当てずっぽうではなく、根拠のある候補から始められるからです。

数値で見た結論

このAIモデル上では、計算箱を243原子の9×9から432原子の12×12へ広げても、CDWへつながる主な「虚振動」は約i82.43 cm-1で変わらず、原子の振付も一致度1.000でした。これはV1が扱う有限距離の相互作用の範囲で、主要結果が箱の大きさに対してほぼ変わらなかったという意味です。DFTや実物の収束を証明するものではありません。

さらにH100上のV1では、箱、振幅、位相、数値対照を含む388構造を直接評価しました。「振動の警報」「エネルギーの下がり方」「その方向へ押す力」という三つの見方は0.08%以内で一致しました。選んだ一本の経路では、負側Q=-0.03303 Aと正側+0.02996 A付近がV1上の谷候補です。続くQuantum ESPRESSOでも、原点近くでは外へ押していた力が、遠い確認点では戻す向きへ反転しました。これは全原子を自由に緩和したCDW構造ではありませんが、「虚数が出た」で止めず、機械学習モデルの予測をDFTで一段深く検査できたことが今回の前進です。

ここで重要なのは、最初からきれいな答えが出たわけではないことです。

• 元の構造には無視できない力が残っていた
• 3×3は有限変位フォノンのforce constantsを作るには小さすぎた
• 6×6はV1内の一部の点では合って見えたが、経路全体の数値検査では不合格だった
• そのため9×9、さらに12×12まで計算した
• 広い調和範囲は不合格となり、狭い範囲で再検証した
• 別の計算手法、Quantum ESPRESSOでもQ=0、±0.010 Aの局所負曲率を確認した
• さらにV1候補Q*正負をQuantum ESPRESSOで照合し、両branchで力の反転を確認した

この「疑い、計算を増やし、条件を一つずつ落とさない」過程こそ、H100が研究にもたらす価値です。

そもそも、なぜ原子の振動を調べるのか

固体の中の原子は、止まっているように見えても振動しています。その振動には、ばねのように元の位置へ戻る安定な動きと、少し押すと別の構造へ進んでしまう不安定な動きがあります。

後者が計算に現れたものを「虚振動」と呼びます。

たとえるなら、ボールが谷底にあるときは、少し押しても戻ります。しかし山の頂上にあるボールは、わずかなきっかけで転がり落ちます。虚振動は、「今の原子配置は谷底ではなく、ある方向に対して山または鞍点にいる」と教えてくれます。

NbSe2では、電子の濃淡と原子のずれが周期的に並ぶCDWと超伝導が、単層まで共存します。CDWは層数、ひずみ、ゲート、積層、基板に敏感で、競合構造の差はsub-meV/atom級です。だから一度だけ計算して答えを出すより、サイズ、変位幅、位相、計算法を変えた対照試験が重要になります。どの原子がどの方向へ動きやすいかを調べることは、CDW構造が生まれる入口を見つけることです。

なぜそれが重要なのでしょうか。量子材料では、わずかな原子の並び替えが電子の流れ方や、CDWと超伝導の関係に結びつきます。どの動きが構造変化の入口になるかを理解することは、将来、層数、ひずみ、電圧、基板などで材料の状態を制御する研究の土台になります。今回の計算だけでデバイス性能を予測したわけではありませんが、その前段となる「どこを詳しく調べるべきか」を具体化しています。

本稿で示す価値はCDWそのものの発見ではありません。公開V1を独立に再現・数値監査し、H100で検証条件を広げ、計算コストの高いDFTへ渡すまでの証拠を切れ目なく残した点にあります。前稿のMD trajectoryとは異なり、今回は静的な調和フォノンとfrozen-mode single-point走査です。

なぜ検証条件を増やすと、計算が急に重くなるのか

この種の計算が従来の方法で不可能だった、あるいは誰も行ってこなかった、という意味ではありません。DFTを使えば、一つ一つの候補を電子の状態から高精度に調べられます。ただし、丁寧に疑うほど計算数が急増することが障壁になります。

小さい計算箱で一度だけ調べても、もっともらしい数字は出ます。しかしNbSe2で競う構造のエネルギー差は非常に小さく、次の条件が結果を変える可能性があります。

• 計算箱をどこまで大きくするか
• 原子をどの方向へ、どれだけ動かすか
• 波を格子のどこから始めるか
• 正負両方向を見るか
• 別の計算法でも同じ傾向になるか

条件を一つ増やすたびに、必要な原子配置も増えます。しかも「きれいな結果が出た点」だけでなく、不合格になった条件や比較用の対照も計算しなければ、偶然の一致を見抜けません。今回の27原子DFT single-pointでさえ、64 MPI ranksで1点約19～38分かかりました。さらに大きな構造や多数の変位を同じ精度で調べれば、計算予算は一層大きくなります。

一方、機械学習モデルだけでは小さな谷を確定できません。そこで今回は、高速なモデルで検証範囲を広げ、精密なDFTを結論に効く少数点へ集中する二段階に分けました。この二段階の狙いは「DFTを使わない」ことではなく、DFTを使う場所を、より多くの対照計算から選べることです。

小さい計算で数字が出ても、すぐには信じない

最初の3×3計算でも虚振動は見えました。しかし、そのまま発表には使いませんでした。

第一の理由は、元の高対称構造で各原子に最大0.1788 eV/Aの力が残っていたことです。全原子の力を足すとゼロでも、上下の原子が逆方向に強く押されていれば、構造は平衡ではありません。

そこで結晶の対称性を保ちながらSe原子の高さを調整しました。SeとNbの高さは1.677975 Aから1.665267 Aへ変わり、V1のforce component最大値は1.87e-8 eV/Aまで下がりました（原子force vector最大値は2.33e-8 eV/A）。

第二の理由は、計算箱の大きさです。原子間の影響を10 A先まで扱うモデルに対して、小さい箱では周期的にコピーされた原子の影響が混ざります。そこで6×6、9×9、12×12へ広げました。

6×6と9×9は、調べた代表点だけならよく一致しました。ところがGamma–M–K–Gammaという経路全体を比べると、最大14.818 cm-1の差と符号の違いが見つかりました。6×6だけで止めていたら、見逃した問題です。

9×9と12×12では、V1上の安定な枝の最大差が0.153 cm-1、不安定な枝でも0.493 cm-1でした。主な虚振動の周波数と原子の動く向きは一致しました。このV1内部の数値検査から、今回のV1調和計算で使う標準条件を9×9、原子を動かす有限変位幅を0.005 Aとしました。ただし元論文では、長距離force constantsを必要とする完全なフォノン収束には大cell学習dataがさらに必要だと説明されています。

またV1上ではM点に弱い虚周波数も現れますが、これは主CDWモードと同格には扱いません。元論文では、より大きいDFT cellでM点が正へ変わることが示されており、DFT確認前のM点はMLIP由来のartifact候補です。

機械学習モデルの数字を、別の計算でも見る

今回のV1は、第一原理計算データを学習した振動特性向けの機械学習ポテンシャルです。前稿でMDに使ったM1とは別の重み・別の用途です。高速ですが、機械学習モデルだけで話を閉じるべきではありません。

まず、モデルが学習に使わず残していた92構造で、いわば「模擬試験」を行いました。

• エネルギー平均誤差：0.948 meV/atom
• 力の成分平均誤差：0.01487 eV/A

これはモデル内部の未使用データに対する結果です。模擬試験に合格しても、未知の構造すべてに正しく答えられる保証にはなりません。

今回V1が予測したCDW方向のenergy低下は1原子あたり約0.52–0.54 meVで、validationのenergy平均誤差0.948 meV/atomより小さい値です。平均誤差をそのまま個々の相対energyの誤差棒にはできませんが、V1だけで谷の深さや正負branchの順位を確定してはいけない、という重要な警告になります。

次に、V1で緩和した1×1構造をQuantum ESPRESSO 7.5で計算しました。64 CPUコアを使い、vdW-DF2-c09、27×27×1 k点、70/560 Ryという条件です。最大残留力は0.007474 eV/Aで、計算前に決めた0.02 eV/Aのscreening gateを満たしました。

つまり、高対称構造の「傾きが大きく残ってはいない」というscreeningには通りました。ただしDFT平衡構造や高精度DFT停留点の証明ではなく、Se高さのDFT緩和が別に必要です。

さらに、虚振動の方向にQ=0、±0.010 Aだけ動かした27原子構造を、同じfunctional・cutoff・smearingと等価なk点密度で計算しました。結果は事前登録した全条件でPASSです。

• Q=-0.010 A：高対称点より3.652 meV/3×3セル低下
• Q=+0.010 A：高対称点より3.752 meV/3×3セル低下
• エネルギーから求めた曲率：-8.226 eV/(式単位·A²)
• 力から求めた曲率：-7.496 eV/(式単位·A²)

両側でエネルギーが下がり、力も原点から外向きでした。つまりV1内部だけでなく、fresh QE/DFTも選んだsingle-q方向に沿って高対称点が局所的に不安定だと判定しました。別実装でraw出力から全27原子の力と曲率を再計算しても、数値と合否は一致しました。

CPU 64 MPI ranksでのQE wall timeは、Q=0が19分02秒、負側が37分36秒、正側が36分28秒でした。H100上のV1による388構造・50秒とは精度階層も計算法も違うため、これを速度倍率にはしません。意味があるのは、GPUで広く探してDFTを判別力のある少数点へ集中できたことです。

山を下見してから、候補地点を地上で測る

H100上のV1探索は、広い山を上空から見渡し、尾根や谷らしい場所を素早く探す「下見」に似ています。Quantum ESPRESSOのDFT計算は、候補地点へ測量班を送り、地上の精密機器で高さと坂の向きを測り直す作業です。

上空からの下見だけでは、本当の谷底だと証明できません。反対に、最初から地上測量だけで多数の候補と対照条件を回ると、時間と計算資源を大きく使います。今回はH100上のV1で対照を含む388構造を調べ、その情報から判別力のある5点をDFTで照合しました。

したがって、ここでのインパクトは「DFTがH100より何倍遅い」という速度競争ではありません。広い探索と精密な確認を組み合わせ、限られたDFT計算を結論に効く場所へ集中できたことです。

「虚振動」を、実際にその方向へ押してみた

虚振動が見つかっただけでは、「少し動かせば本当にエネルギーが下がるのか」「どこで止まるのか」は分かりません。そこで、q_CDWの原子の動きに沿って、振幅と位相を変えた構造を実際に作り、エネルギーと力を直接計算しました。

平たく言えば、q_CDWは原子集団の「振付と繰り返し周期」、Qはその振付で原子をどれだけ大きく動かすか、phaseは格子のどこから振付を始めるかです。厳密には、q_CDW=(1/3,-1/3,0)はprimitive reciprocal basisで2/3 Γ–Mに相当します。Qは一般的な質量重み付きnormal coordinateではなく、原子1個当たりCartesian変位vectorのRMS（単位Å）です。曲率はこの規格化に依存します。式単位（f.u.）はNbSe2の3原子です。phase 0はunit-cell originと固有vectorのglobal phaseを固定した数値規約で、絶対的な観測位相ではありません。

最初の広域scanは、3×3、6×6、9×9、12×12、0を含む9段階の振幅（非零は正負）、6位相を組み合わせた、対照を含む388構造評価です。H100×2のSlurm実測は50秒でした。

ここでも、最初から合格にはしませんでした。0.005、0.010、0.020 Aを一つの「小さな変位」として扱う事前基準を決めていましたが、0.020 Aではすでに曲線が単純な放物線から外れ始めていました。そのため広域scanの正式判定はFAILのまま残しています。

範囲を0.0025、0.005、0.010 Aへ狭めると、三つの方法が同じ答えを出しました。

• 虚周波数から計算した曲率：-6.64496 eV/(式単位·A²)
• energyの下がり方から計算した曲率：-6.63992
• forceの向きから計算した曲率：-6.64386

差は最大でも0.08%です。つまり「振動数が虚数になった」「実際にenergyが下がった」「forceがその方向へ押した」という三つの見方が、同じV1モデルの中で一本につながりました。

さらに細かく調べると、phase 0の一次元経路では、負側Q=-0.03303 Aと正側Q=+0.02996 A付近で、その方向の傾きがほぼゼロになりました。各点の左右±0.001 Aも新しく計算し、中心のenergyが両隣より低いこと、左右のforceの向きが反転すること、曲率が正になることを確認しました。

高対称構造からのenergy低下は、負側1.609 meV/式単位、正側1.561 meV/式単位です。正負で少し違うのは計算誤差とは限りません。この波は3回繰り返すと結晶周期に戻るため、三次の「lock-in」が許され、左右対称でなくてもよいからです。

これはV1上の一次元局所極小候補です。固定したcell・一つのq・一つの位相という道筋の上だけの谷であり、道筋と直角な方向にはまだ力が残っています。したがって「CDW構造を発見した」「基底状態を確定した」とは書きません。

DFTでも谷を挟めたか

V1で得た二つの候補を、Stage Aと同一のQuantum ESPRESSO条件でsingle-point計算しました。Q=0と±0.010 Aを含めると、DFTに使ったのは合計5点です。ここで比較するのは、QE内でQ=0を基準にしたrelative free energyだけです。V1とQEのabsolute energyは混合しません。

結果を先に平たく言うと、谷の内側ではボールを外へ押す坂、少し先ではボールを戻す坂になっていました。坂の向きが二地点の間で逆になるため、その間に傾きがゼロになる場所が少なくとも一つあります。この「二地点で谷の範囲を挟んだ」という判定を、以下ではBRACKETEDと呼びます。

QE Stage Bの正式statusはPASS、decision_codeはPASS_DFT_LOWERED_FROZEN_MODE_SEEDです。

負側branchのforce分類はBRACKETED、正側branchのforce分類はBRACKETEDです。

負側Q*のQE相対自由energyは-17.123 meV/3×3 cell、正側Q*は-14.244 meV/3×3 cellです。1式単位当たりでは、それぞれ-1.903、-1.583 meVです。どちらも事前基準の-1 meV/3×3 cellを大きく超えてQ=0より低くなりました。

力はどうでしょうか。mode方向のgeneralized forceをG=(1/9)ΣF_i·p_i、エネルギーの傾きをD=dF/(f.u.)/dQ=-Gと定義します。small-Qでは原点から外向きだったDが、負側Q*で-0.01235、正側Q*で+0.02348 eV/(Å·f.u.)となり、両側とも復元向きへ反転しました。Q=0の小さな残差を差し引いても結論は同じです。

これは、負側ではQ=-0.010から-0.03303 Aの間、正側ではQ=+0.010から+0.02996 Aの間に、傾きがゼロとなる一次元stationary regionが少なくとも一つある、という意味です。Q*そのものが停留点だとは限りません。イメージとしては、山頂から転がり出したボールが、調べた先では「戻る向き」に押されていたため、その途中に谷底方向の折り返しを挟めた、ということです。

Stage Bは固定cellのfrozen-mode照合です。single-pointをDFT極小とは扱いません。原子を全自由度で緩和した局所極小、競合CDW構造との順位、global minimumを示すには、別にDFT構造緩和が必要です。

H100×2が変えたのは、計算時間より「疑える回数」

今回の成果は、CDWを新発見したことでも、388点すべてをDFTで証明したことでもありません。3×3を疑い、6×6を不合格にし、9×9と12×12の一致を確認し、広すぎた調和範囲も不合格のまま記録したうえで、V1の予測を別枠の27原子DFT 5点へつないだことです。高速化によって増えたのは、都合のよい答えではなく、誤りを落とすための検査回数でした。

2枚のH100は別々の候補構造を同時に担当しました。一つの構造を2枚で半分の時間にした測定ではありません。対照を含む388構造を50秒で評価できたため、不合格になった範囲を見た後に、範囲を狭めた再計算と局所三点検査を別runで追加できました。

研究では一回の最速記録より、「箱が小さいせいではないか」「原子の動かし幅のせいではないか」を追加計算で疑えることが結論の強さにつながります。H100×2の価値は、結果を疑い、次の条件を考える機会を増やすことです。

まだ言えないこと

今回確認したのは、検証したV1ポテンシャル面上で0 Kの高対称単層NbSe2に強い調和不安定性があり、fresh QE/DFTでも選んだsingle-q座標の原点が局所的に不安定で、正負両branchに一次元stationary regionを挟めたことです。

まだ、次のことは確定していません。

• 最終的にどのCDW構造が最も安定か
• 全原子自由度を緩和したとき、どの局所極小へ到達するか
• CDW転移温度が何Kか
• 超伝導転移やその機構

調和フォノンや短い古典MDだけではCDW転移温度T_CDWを定量化できません。次段階は、同じDFT referenceに対するsingle-q正負と既知のhollow/filled/hexagonal・3q位相seedの全自由度緩和です。その後、非調和性と核量子効果を扱うV1+SSCHAへ進み、SSCHA ensembleから代表構造をDFT spot-checkし、一次転移の可能性があればphase間free-energyも比較します。SSCHAは唯一の方法ではなく、超伝導転移温度T_cを直接求める計算でもありません。T_cにはDFPT/EPWで電子–フォノン結合とα²Fを求め、Coulomb・SOC等の仮定を置いた別工程が必要です。

お客様の計算でも、「どこまで疑えるか」を測ります

GPU導入前に知りたいのは、単純な最高速度だけではないはずです。

• 今のCPU計算は何時間かかっているか
• 同一model・code・workloadでCPUからGPUへ移したとき、同じ入力を何本比較できるか
• 1枚と2枚で、独立した条件をどう分担できるか
• 同一workloadの実測で、速度だけでなく精度確認や収束試験まで含めて何日短くできるか

ServerGearでは、お客様の入力を使うおまかせ性能測定をご用意しています。H100単体とH100×2搭載システムの相談は、hanabusa@server-gear.comまでお寄せください。

測定上の注意

• H100はPCIe 80 GB×2、FP64、TF32 offです。
• 調和フォノンはV1+ASE+Phonopy 4.4.0。QE確認は別ホストのCPU 64 MPI ranksです。
• LAMMPSは代表構造1点の独立実装検査に使用し、Python/ASE経路とのenergy差は0、force component最大差は7.83e-14 eV/Aでした。今回の388構造探索とフォノン計算の主経路ではなく、MDも実行していません。
• QEの約35秒は3原子の1×1残留力screeningで、H100側は最大432原子の有限変位計算です。系の大きさも目的も異なるため、両時間を割った速度倍率は示していません。
• QE Stage Aの27原子single-pointはQ=0が19分02秒、負側37分36秒、正側36分28秒です。これもH100上のV1とのhardware speedupには使いません。
• QE Stage Bの27原子single-pointは負側Q*が38分00秒、正側Q*が37分20秒です。これらは候補点のDFT照合時間で、GPUとの速度倍率にはしません。
• q_CDW広域scanの対照を含む388構造評価・50秒は、V1によるGPU force evaluationです。DFT 388点を50秒で計算したという意味ではありません。
• 2枚のH100はconfigurationを分担するtask parallelismです。一つのconfigurationを2 GPUへ分割したstrong-scaling値ではありません。
• M点の弱い虚周波数はDFT未確認のMLIP artifact候補であり、物理結論には使っていません。
• 0.0025/0.005 Aの主モードは収束しましたが、ほぼゼロの弱い枝の交差位置には最大約1.99 cm-1の変位幅依存性があります。
• 元論文：npj Computational Materials (2026)