Machine learned features from density of states for accurate adsorption energy prediction

# ref
https://www.nature.com/articles/s41467-020-20342-6

# 概要
畳み込みニューラル ネットワーク ベースの機械学習モデルを使用してエネルギーを予測し、電子状態密度 (DOS) から主要な特徴を自動的に取得するアプローチを示す。この研究では、吸着エネルギーの予測に焦点を当てる。

### 吸着エネルギーとは...?
吸着エネルギーとは、物質が固体表面に付着する際に関与するエネルギーのことを指す。吸着エネルギーは、物質が固体表面にどの程度強く結合するかを示す指標であり、接着現象や濡れ性、粒子分散性に影響を与える。

計算方法は、「表面・分子間の相互作用エネルギーをCounterpoint法で補正したエネルギー」と、「表面・分子の吸着状態と孤立状態のエネルギー差」からエネルギーを求める。

# 手法
用いる機械学習モデルはCNN。DFT(密度関数理論)から得られた特性をCNNに供給し、吸着エネルギーを予測するための新しい特徴量を抽出。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3487085/1d46b71f-0912-9799-8a27-3598cb787b7c.png)


# 実装
### 機械学習モデル
CNN

### CNNモデル構成

#### 畳み込み層
以下はデフォルトの値
```python
def dos_featurizer(channels):
    input_dos = Input(shape=(2000, channels))
    x1 = AveragePooling1D(pool_size=4, strides=4, padding="same")(input_dos)
    x2 = AveragePooling1D(pool_size=25, strides=4, padding="same")(input_dos)
    x3 = AveragePooling1D(pool_size=200, strides=4, padding="same")(input_dos)
    x = Concatenate(axis=-1)([x1, x2, x3])
    x = Conv1D(50, 20, activation="relu", padding="same", strides=2)(x)
    x = BatchNormalization()(x)
    x = Conv1D(75, 3, activation="relu", padding="same", strides=2)(x)
    x = AveragePooling1D(pool_size=3, strides=2, padding="same")(x)
    x = Conv1D(100, 3, activation="relu", padding="same", strides=2)(x)
    x = AveragePooling1D(pool_size=3, strides=2, padding="same")(x)
    x = Conv1D(125, 3, activation="relu", padding="same", strides=2)(x)
    x = AveragePooling1D(pool_size=3, strides=2, padding="same")(x)
    x = Conv1D(150, 3, activation="relu", padding="same", strides=1)(x)
    shared_model = Model(input_dos, x)
    return shared_model
```
- 入力として受け取るのは形状が(2000,channels)のDOSデータ。
- 結果的な出力の形状は(4, 150)となる。
1. AveragePooling1D（x1, x2, x3）
異なるpool_sizeで平均プーリングされる。strides=4なので長さは1/4に。
    > input:(2000,channels)
    output:x1(500, channels), x2(500, channels), x3(500, channels)
1. Concatenate
    x1, x2, x3を結合している。
    > input:x1(500, channels), x2(500, channels), x3(500, channels)
    output:(500, 3*channels)
1. Conv1D（50フィルタ、20カーネルサイズ）
    strides=2なので長さが半分に。フィルタ数が50なのでチャンネル数は50に。
    学習可能パラメータ数は50*(20*3channels+1)
    > input:(500, 3*channels)
    output:(250, 50)
1. BatchNormalization
    データ形状は変化なし
    各入力チャネルに対してスケールパラメータとシフトパラメータがあるため、学習可能パラメータは2(50+75+100+125+150)←2*入力チャネル数
    > input:(250, 50)
    output:(250, 50)
1. Conv1D（75フィルタ、3カーネルサイズ）
    strides=2なので長さが半分に。フィルタ数が75なのでチャンネル数は75に。
    学習可能パラメータ数は75*(3*50+1)
     > input:(250, 50)
    output:(125, 75)
1. AveragePooling1D
    strides=2なので長さが半分に。
    > input:(125, 75)
    output:(63, 75)
1. Conv1D（100フィルタ、3カーネルサイズ）
    strides=2なので長さが半分に。フィルタ数が100なのでチャンネル数は100に。
     学習可能パラメータ数は100*(3*75+1)
    > input:(63, 75)
    output:(32, 100)
1. AveragePooling1D
    strides=2なので長さが半分に。
    > input:(32, 100)
    output:(16, 100)
1. Conv1D（125フィルタ、3カーネルサイズ）
    strides=2なので長さが半分に。フィルタ数が125なのでチャンネル数は125に。
    学習可能パラメータ数は125*(3*100+1)
    > input:(16, 100)
    output:(8, 125)
1. AveragePooling1D
    strides=2なので長さが半分に。
    > input:(8, 125)
    output:(4, 125)
1. Conv1D（150フィルタ、3カーネルサイズ）
    strides=1なので長さは変化なし。フィルタ数が150なのでチャンネル数は150に。
    学習可能パラメータ数は150*(3*125+1)
    > input:(4, 125)
    output:(4, 150)

#### モデル
```python
def create_model(shared_conv, channels):

    ###Each input represents one out of three possible bonding atoms
    input1 = Input(shape=(2000, channels))
    input2 = Input(shape=(2000, channels))
    input3 = Input(shape=(2000, channels))

    conv1 = shared_conv(input1)
    conv2 = shared_conv(input2)
    conv3 = shared_conv(input3)

    convmerge = Concatenate(axis=-1)([conv1, conv2, conv3])
    convmerge = Flatten()(convmerge)
    convmerge = Dropout(0.2)(convmerge)
    convmerge = Dense(200, activation="linear")(convmerge)
    convmerge = Dense(1000, activation="relu")(convmerge)
    convmerge = Dense(1000, activation="relu")(convmerge)

    out = Dense(1, activation="linear")(convmerge)
    # shared_conv.summary()
    model = Model([input1, input2, input3],out)
    return model
```
1. ３つの金属原子のDOSをinputし、さっきのdos_featurizer関数で特徴量抽出。`(2000,channels)→(4,150)`
1. それらを結合し、平坦化`(4,150)*3→(1800,1)`
1. ２０％をドロップアウト。`(1800,1)`
1. 全結合層4つ
    `(1800,1)→(200,1)→(1000,1)→(1000,1)→(1,1)`
1. 答えと比較し、誤差を出す。
1. 逆伝播し、全結合層の各ニューロンの重みを更新
1. 繰り返し
1. 各エポックごとにvalidataion_dataを出す。(現時点では出してるだけ)

その後、`model.predict`で予測が行われる。予測結果は形状が`(サンプル数, 1)`の2次元配列となり、これを`.reshape(len(predicted_data))`で1次元配列に変換。したがって、最終的な予測結果`train_out, test_out`は形状がそれぞれ`(サンプル数,)`の1次元配列となる。

入力から予測までの流れは下の画像の通り。
![Screenshot from 2023-11-06 10-25-53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3363422/2864cfa9-1738-fb3f-9eb7-475e0a1cbe40.png)


### データ
#### 入力
金属とそれに吸着させる原子のDOS。基本的には金属のDOSを原子３つ分。`malti_absorbut`のオプションで、吸着物質のデータをinputするかどうかを設定できる。

これらのDOSは、各々が異なるチャンネルとして扱われ、それぞれが独立した情報を提供。最大3段階の原子を入力に含めることができる。原子一つに付き９チャネルとなり、原子が増えるごとに18,27とチャネル数が増えていく。各DOSデータは2000の長さを持つ1次元配列。

形状は(サンプル数, DOSファイルの長さ, チャネル数)のnumpy配列として保存される。
平均をゼロ、すべてのチャネルにわたり分散を１にスケーリングして標準化。

37000の吸着エネルギーデータ。
合金は、化学量論比0,0.25,0.5,0.75,1の37種類の遷移金属元素と非繊維金属元素で構成されており、2000種類ある。これらの合金表面に対し、様々な吸着物質、吸着位置、表面構造などがあり、それらの総数が37000となる。なお、吸着物質はH, C, N, O, Sとそれらの水酸化化合物である。(具体的にはCH、CH₂、CH₃、NH、OH、SH )

#### 出力
予測した吸着エネルギーと実際の吸着エネルギーの差(eV)


#### 評価指標
予測した吸着エネルギーと実際の吸着エネルギーの絶対誤差(eV)


# 再現

1. 下記のgitからコードをclone
`git clone https://github.com/vxfung/DOSnet`
2. 必要パッケージのインストール
`pip install -r requirements.txt`
1. 以下urlからzipファイルをダウンロード&解凍(解凍場所はDOSnet/)
https://figshare.com/articles/dataset/DOSnet_data/14511978
`tar -xvf "your_file_here.tar.gz"`
1. Main.py line156を下記に書き換え
`model = Model([input1, input2, input3],out)`
1. 実行
`python Main.py --multi_adsorbate=0 --data_dir='CH_data' --save_model=1 --batch_size=128 --channels=9`
1. 出力
`python3 plot.py --dataname="CH_data"`
結果
train MAE:  0.15028492
train RMSE:  0.19413818188977927
test MAE:  0.19662598
test RMSE:  0.2568728779382601

CHデータをDOSnetに通した結果をプロット
散布図において、横軸が実際の値で、縦軸が予測値である。
また、散布図外のヒストグラムはそれぞれの一定範囲ごとの予測値の数を表している。
![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3363422/489d5225-a470-4178-8d99-84b21866c5f9.png)