はじめに

前回の記事では、MODNet (Material Optimal Descriptor Network) とは何か、そしてその理論と応用例について解説しました。 MODNetは、小規模なデータセットでも高精度な材料特性予測を可能にする優れたモデルです。

今回はMODNetを実際にインストールし、簡単なサンプルコードを実装して材料特性予測を行うための手順を解説します。 MODNetはオープンソースのPythonパッケージとして公開されており、誰でも簡単に利用することができます。

実行環境

筆者の実行環境は以下になります。今回はanacondaの利用を前提とします。

Windows 11
conda 24.11.3
jupyter 1.1.1

インストール

公式のGithubに従って環境構築とインストールを行っていきます。

Githubからパッケージをダウンロードする

git cloneを使ってリポジトリからパッケージをダウンロードします。

コマンドプロンプト

git clone git@github.com:ppdebreuck/modnet

ダウンロードしたmodnetにディレクトリを移動します。

コマンドプロンプト

cd modnet

仮想環境の構築

anacondaで仮想環境を作ります。pythonのバージョンは3.9を指定します。

コマンドプロンプト

conda create -n modnet python=3.9

そして、作った仮想環境をアクティベートします。

コマンドプロンプト

conda activate modnet

pipによるインストール

pipを使ってパッケージのインストールを行います。requirements.txtを使って依存パッケージをインストールし、editableオプションを使ってmodnetをインストールします。

コマンドプロンプト

pip install -r requirements.txt
pip install -e .

これでインストールは完了です。

MODNetの実装

それでは、jupyterを使ってMODNetの実装を行っていきます。今回は材料の組成を入力としてバンドギャップを予測するというタスクを想定します。こちらもGithubにガイドがあるため、それを参考にして実装します。

必要モジュールのインポート

まず最初に必要となるモジュールのインポートを行っていきます。matminerはデータセットのダウンロードや材料の特徴量生成をサポートするモジュールです。今回は詳細について触れませんが、別で記事を書くかもしれません。

python

import numpy as np
from matminer.datasets import load_dataset
from modnet.models import MODNetModel
from modnet.preprocessing import MODData
import matplotlib.pyplot as plt 
from pymatgen.core import Composition
from matminer.datasets import load_dataset

データセットのダウンロード

matminerのload_dataset関数を使用してデータセットのダウンロードを行います。そして、後の特徴量生成のため、組成をpymatgenのComposition型に変更します。

python

df = load_dataset("matbench_expt_gap")
df["composition"] = df["composition"].map(Composition) # Compositionオブジェクトへマッピングする
df.head() # データフレームを確認

実行結果は次のようになります。材料の組成と、対応するバンドギャップが得られていることが分かります。

ここで、バンドギャップがどのように分布しているかを確認してみましょう。matploblibを使ってプロットします。
このデータセットはバンドギャップがゼロの材料も含まれるため、それを分けて表示します。

python

fig, ax = plt.subplots(facecolor="w")
ax.hist(df.where(df["gap expt"] == 0)["gap expt"], bins=1, density=False, label="Zero band gap")
ax.hist(df.where(df["gap expt"] > 0)["gap expt"], bins=25, density=False, label="Non-zero band gap")
ax.set_ylabel("Frequency")
ax.set_xlabel("Band gap (eV)")
ax.legend()

実行結果は以下のようになります。オレンジが非ゼロのバンドギャップの分布です。2.0eVあたりをピークに、右に裾を引くような分布になっているようですね。

MODDataオブジェクトの作成

MODNetは、主にMODDataクラスとMODNetModelクラスによって構成されます。MODDataクラスの方は学習に使用するデータセットを操作するためのクラスです。最初にダウンロードしたデータセットは組成のみの情報でしたが、このクラスに渡すことで様々な情報を取得して特徴量生成などを行うことができます。引数であるmaterialsにはcompositionのカラム、そしてtargetsには目的変数となるバンドギャップのカラムを渡します。また、target_namesでは目的変数に名前を付けることができ、これは自由に変えられます。

python

data = MODData(
    materials=df["composition"], # 作成したデータフレームのcommpositionカラムを渡す
    targets=df["gap expt"], 
    target_names=["gap_expt_eV"]
)

特徴量生成

MODDataオブジェクトを使って特徴量生成を行います。featurize関数を呼び出すだけでOKです。

python

data.featurize()

すると、このように特徴量生成が進行します。環境にもよるかと思いますが結構時間がかかります。筆者の環境では10分ほどかかりました。ここでは特徴量選択はまだ行わず、生成できる限りの特徴量が計算されるようです。

実行結果

2025-03-11 22:18:46,442 - modnet - INFO - Computing features, this can take time...
2025-03-11 22:18:46,445 - modnet - INFO - Applying composition featurizers...
2025-03-11 22:18:46,457 - modnet - INFO - Applying featurizers (AtomicOrbitals(), AtomicPackingEfficiency(), BandCenter(), ElementFraction(), ElementProperty(data_source=<matminer.utils.data.MagpieData object at 0x0000015022757AF0>,
                features=['Number', 'MendeleevNumber', 'AtomicWeight',
                          'MeltingT', 'Column', 'Row', 'CovalentRadius',
                          'Electronegativity', 'NsValence', 'NpValence',
                          'NdValence', 'NfValence', 'NValence', 'NsUnfilled',
                          'NpUnfilled', 'NdUnfilled', 'NfUnfilled', 'NUnfilled',
                          'GSvolume_pa', 'GSbandgap', 'GSmagmom',
                          'SpaceGroupNumber'],
                stats=['minimum', 'maximum', 'range', 'mean', 'avg_dev',
                       'mode']), IonProperty(), Miedema(ss_types=['min'], struct_types=['inter', 'amor', 'ss']), Stoichiometry(), TMetalFraction(), ValenceOrbital(), YangSolidSolution()) to column 'composition'.
MultipleFeaturizer: 100%
 4604/4604 [06:19<00:00,  8.21it/s]
2025-03-11 22:25:13,216 - modnet - INFO - Data has successfully been featurized!

このfeaturize関数について、中でどのような処理が動いているのか少し見てみましょう。
まず、MODDataクラスのオブジェクトを作成した際に、自動で入力データに応じたFeaturizer(特徴量生成器)が選択されるようになっているようです。クラスのコンストラクタを見てみると、今回はこのelifの部分に入り、self.featurizerにFEATURIZER_PRESETS[DEFAULT_COMPOSITION_ONLY_FEATURIZER]が代入されています。

これらの変数はどちらも別ファイルで定義されているもので、FEATURIZER_PRESETSは辞書型のオブジェクトで、DEFAULT_COMPOSITION_ONLY_FEATURIZERは文字列として定義されており、この文字列によりCompositionOnlyMatminer2023Featurizerというfeaturizerが選択される仕組みになっています。

MODDataクラスからfeaturize関数を呼び出すと、最終的にはここで参照されたCompositionOnlyMatminer2023Featurizerクラスの親クラスであるMODFeaturizerというクラスのfeaturize関数が呼び出され、最終的には_fit_apply_featurizers関数が呼び出されます。そして、この関数の中ではMatminerで用意されたクラスであるMultipleFeaturizerが使用され、featurize_dataframeという関数によって最終的な特徴量生成が行われます。

このようにかなり複雑なフローになっているので追いかけるのが大変でしたが、基本的には自動的に入力のデータ構造を読み込み、Matminerを使って適切な特徴量を生成してくれるようなので、ユーザーレベルではここまで知る必要はないかもしれません。

この処理により生成された特徴量を見てみます。df_featurizedという属性に特徴量がデータフレームとして入るので、このデータフレームの形とカラムを表示します。

python

print(data.df_featurized.shape)
print(data.df_featurized.columns)

実行結果

(4604, 270)
Index(['AtomicOrbitals|HOMO_character', 'AtomicOrbitals|HOMO_element',
       'AtomicOrbitals|HOMO_energy', 'AtomicOrbitals|LUMO_character',
       'AtomicOrbitals|LUMO_element', 'AtomicOrbitals|LUMO_energy',
       'AtomicOrbitals|gap_AO',
       'AtomicPackingEfficiency|mean simul. packing efficiency',
       'AtomicPackingEfficiency|mean abs simul. packing efficiency',
       'AtomicPackingEfficiency|dist from 1 clusters |APE| < 0.010',
       ...
       'ValenceOrbital|avg s valence electrons',
       'ValenceOrbital|avg p valence electrons',
       'ValenceOrbital|avg d valence electrons',
       'ValenceOrbital|avg f valence electrons',
       'ValenceOrbital|frac s valence electrons',
       'ValenceOrbital|frac p valence electrons',
       'ValenceOrbital|frac d valence electrons',
       'ValenceOrbital|frac f valence electrons',
       'YangSolidSolution|Yang omega', 'YangSolidSolution|Yang delta'],
      dtype='object', length=270)

このように、270個の特徴量が生成されていることが分かります。名前から推測する限りでは電子状態に関するものが多く生成されているようですね。

データセットの分割

次に、学習を行うためにデータセットの分割を行います。scikit-learnのtrain_test_split関数を用いて行います。

python

from sklearn.model_selection import train_test_split
split = train_test_split(range(100), test_size=0.1, random_state=1234)
train, test = data.split(split)

特徴量の選択

前回の記事で書いたように、MODNetの大きな利点は特徴量選択にあり、これにより少量のデータセットでも効果的に学習を行うことができます。特徴量の選択はfeature_selectionという関数を使用します。今回は270個の特徴量を100個に絞ってみます。

python

train.feature_selection(n=100) #100個の特徴量を選択する

ソースコードを見てみると、preprocessing.pyの get_features_dynという関数の中で正規化相互情報量(NMI)やRRスコアの計算が行われているようです。

選択された特徴量はoptimal_features_by_targetという属性に入ります。これは辞書型のオブジェクトであり、MODDataオブジェクトを作成するときに指定した"gap_expt_eV"をキーとして使うことで特徴量の一覧をリストとして取得できます。

python

train.optimal_features_by_target["gap_expt_eV"]

実行結果

['ElementProperty|MagpieData mode NUnfilled',
 'ElementProperty|MagpieData avg_dev NValence',
 'ElementProperty|MagpieData mean Number',
 'ElementProperty|MagpieData range Row',
 'ElementProperty|MagpieData avg_dev MeltingT',
 'AtomicPackingEfficiency|dist from 5 clusters |APE| < 0.010',
 'ElementProperty|MagpieData maximum NdUnfilled',
 'ElementProperty|MagpieData avg_dev Electronegativity',
 'ValenceOrbital|frac p valence electrons',
 'ElementProperty|MagpieData avg_dev NUnfilled',
 'ElementProperty|MagpieData minimum Number',
 'Stoichiometry|2-norm',
 'ElementProperty|MagpieData avg_dev GSvolume_pa',
 'ElementProperty|MagpieData mean MeltingT',
 'ElementProperty|MagpieData mean GSvolume_pa',
 'ElementProperty|MagpieData avg_dev NpUnfilled',
 'ElementProperty|MagpieData avg_dev MendeleevNumber',
 ...以下略

モデルの作成

次に、いよいよモデルの作成を行います。モデルはMODNetModelクラスを用いて作成します。

python

model = MODNetModel([[['gap_expt_eV']]],
                    weights={'gap_expt_eV':1},
                    num_neurons = [[256], [128], [16], [16]],
                    n_feat = 100,
                    act =  "elu"
                   )

ここで、各引数の意味は以下の通りです。

[[['gap_expt_eV']]]: ターゲットとなる変数のリスト。複数の特徴量を選択できるようにネスト化したリストになっている
weights: ターゲートとなる変数に付与する重み。今回は予測する変数が一つなので1ですが、複数の物性を予測する際には重みが大きいほどlossへの寄与が大きくなります。
num_neurons: ニューラルネットの各レイヤーのノード数
n_feat: 入力する特徴量の数。
act: ネットワークの活性化関数。eluはExponential Linear Unitを表します。

このようにモデルのビルド時もn_featで特徴量の数を指定するようになっていますが、実際はモデル学習時に使用するfit関数内で dataのget_optimal_descriptors関数が呼び出され、その時に選定された特徴量から上位n_feat個の特徴量を使うようになっています。なお、n_featより選定した特徴量が少ないとエラーを出力するようです。

モデルの学習

以下のコードでモデルの学習を行います。バリデーションデータの割合、学習率、バッチサイズ、損失関数、エポック数などを指定します。

python

model.fit(train,
          val_fraction = 0.1,
          lr = 0.0002,
          batch_size = 64,
          loss = 'mae',
          epochs = 200,
          verbose = 1,
         )

実行するとこのように学習が進行します。今回は詳しく見ませんが、ある程度学習が収束していることが読み取れます。

実行結果

epoch 0: loss: 1.019, val_loss:0.918 val_mae:0.918
epoch 1: loss: 0.851, val_loss:0.700 val_mae:0.700
epoch 2: loss: 0.746, val_loss:0.570 val_mae:0.570
epoch 3: loss: 0.681, val_loss:0.481 val_mae:0.481
epoch 4: loss: 0.632, val_loss:0.450 val_mae:0.450
epoch 5: loss: 0.616, val_loss:0.431 val_mae:0.431
epoch 6: loss: 0.603, val_loss:0.416 val_mae:0.416
---中略---
epoch 193: loss: 0.154, val_loss:0.431 val_mae:0.431
epoch 194: loss: 0.150, val_loss:0.431 val_mae:0.431
epoch 195: loss: 0.140, val_loss:0.419 val_mae:0.419
epoch 196: loss: 0.140, val_loss:0.423 val_mae:0.423
epoch 197: loss: 0.144, val_loss:0.420 val_mae:0.420
epoch 198: loss: 0.141, val_loss:0.420 val_mae:0.420
epoch 199: loss: 0.136, val_loss:0.406 val_mae:0.406

モデルのテスト

トレーニングしたモデルを使ってテストデータのバンドギャップを予測してみます。

python

pred = model.predict(test)
pred.head() #上から5番目まで確認

実行結果

gap_expt_eV
id40	0.029450
id35	3.807828
id81	-0.235226
id61	2.046163
id98	1.723046

そして、テストデータに対する平均絶対誤差(MAE)を計算します。

python

mae_test = np.absolute(pred.values-test.df_targets.values).mean()
print(f'mae: {mae_test}')

実行結果

mae: 0.3859002037048341

テストデータに対する誤差は0.386eV程度となりました。あまり良い結果とは言えませんね。モデルを改善するにはハイパーパラメータの調整や、どのようなデータに対して精度が悪いのかなど、解析していく必要がありそうです。

まとめ

以上、シンプルなタスクではありますがMODNetを実装し、材料物性の予測を行ってみました。ソースコードの依存関係が複雑だったので、自分でいろいろとカスタマイズして使っていくには少し難しそうだというのが正直な感想でした。ただ、強力なモデルであることには間違いないと思いますので、また機会があれば複数物性の予測などにも挑戦してみたいと思います。間違っている記載などあればご指摘いただければと思います。

参考文献