AENETを使ってみる

#はじめに
[aenet](http://ann.atomistic.net/)を使ってポテンシャルエネルギー曲線の学習方法を学習してみます。
学習に用いたデータは[こちら](https://github.com/ikuhamada/aenet-learn)からダウンロードが可能です。

#準備
##aenetのコンパイル
私の場合``${HOME}/AENET/src``にソースコードを置きコンパイルを行います。

```Bash
mkdir -p ~/AENET/src
```

```Bash
cd ~/AENET/src
```

```
wget http://ann.atomistic.net/files/aenet-2.0.3.tar.bz2
```

```
tar jxvf aenet-2.0.3.tar.bz2
```

``aenet-2.0.3/lib``に移動し、``Makefile``を編集して環境に応じて適当なコンパイラとオプションを指定します。例えば

```Makefile
FC = ifort -c
```
```Makefile
FCFLAGS = -O2 -xSSE4.2 -axCOMMON-AVX512,CORE-AVX512,CORE-AVX2,CORE-AVX-I,AVX
```

などを指定し``make``を実行します。

その後``../src``に移動し、ifortとIntel MPIを使用する場合は
```
make -f makefile/Makefile.ifort_intelmpi
```
を実行してaenet本体のコンパイルを行います。コンパイルに成功すると``bin/``ディレクトリに
- ``generate.x-2.0.3-ifort_intelmpi``
- ``predict.x-2.0.3-ifort_intelmpi``
- ``train.x-2.0.3-ifort_intelmpi``
が作成されます。これらへのエイリアスを``generate.x``、``predict.x``、``train.x``という名前で作成しておきます。

##データのダウンロード
GitHub上のデータを適当なディレクトリでダウンロードします。

```
git clone https://github.com/ikuhamada/aenet-learn.git
```

#ポテンシャルの学習

## レナードジョーンズ二量体

- ディレクトリ``LJ/``

ポテンシャルエネルギー学習の最も簡単な例としてレナードジョーンズ二量体を考えてみましょう。この例はCooperら[Cooper et al., npj Comput. Mater. **6**, 54 (2020)]で紹介された結果を再現してみたいと思います。。

論文に記載されたパラメーターを使ってポテンシャルエネルギーを書いてみます。
![lj.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/837941/bceb989e-a958-46fb-591c-dcce81dc1500.png)
対応する構造はダミー原子として炭素を使用して作成し``xsf/``に置いてあります。

先ずは``01-generate/``で以下の入力ファイル（``generate.in``）を使ってトレーニングセットを作成します。

```generate.in
OUTPUT LJ.train

TYPES
1
C 0.0

SETUPS
C C.fingerprint.stp

FILES
8
../xsf/structure01.xsf
../xsf/structure02.xsf
../xsf/structure03.xsf
../xsf/structure04.xsf
../xsf/structure05.xsf
../xsf/structure06.xsf
../xsf/structure07.xsf
../xsf/structure08.xsf
```
入力ファイルに加えて構造のフィンガープリントを記述するファイル（``SETUPS``）を用意する必要があります。この例ではChebyshev基底を使います。

```C.fingerprint.stp
DESCR
N. Artrith and A. Urban, Comput. Mater. Sci. 114 (2016) 135-150.
N. Artrith, A. Urban, and G. Ceder, Phys. Rev. B 96 (2017) 014112.
END DESCR

ATOM C

ENV 1
C

RMIN 0.5D0

BASIS type=Chebyshev
radial_Rc = 4.0 radial_N = 10 angular_Rc = 4.0 angular_N = 0
```
``generate.x``を実行します。

```bash
generate.x train.in > train.out
```
トレーニングセットは``LJ.train``に出力されます。

次に``02-train/``に移動します。``01-generate/``で作成した``LJ.train``をコピーまたはシンボリックを作成します。

```bash
ln -s ../01/LJ.train
```
以下の入力ファイル``train.in``を使って学習を行います。

```train.in
TRAININGSET LJ.train
TESTPERCENT 10
ITERATIONS  1000

TIMING

METHOD
bfgs

NETWORKS
C C.5t-5t.ann 2 5:tanh 5:tanh
```

``train.x``の実行

```bash
train.x train.in > train.out
```

計算が問題無く終わるとポテンシャルファイル``C.5t-5t.ann``が生成されます。
得られたポテンシャルの検証を``03-test/``で行います。以下の入力ファイル（``predict.in``）を用いてエネルギーの計算を行います。ここではトレーニングに使用した構造を用います。

```predict.in
TYPES
1
C

NETWORKS
 C C.5t-5t.ann 2 5:tanh 5:tanh

FILES
8
../xsf/structure01.xsf
../xsf/structure02.xsf
../xsf/structure03.xsf
../xsf/structure04.xsf
../xsf/structure05.xsf
../xsf/structure06.xsf
../xsf/structure07.xsf
../xsf/structure08.xsf
```
``predict.x``の実行

```bash
predict.x predict.in > predict.out
```
ニューラルネットワークポテンシャルを用いて予測されたエネルギーは``predict.out``に出力されます。``predict.out``のメッセージ

```
 ----------------------------------------------------------------------
                           Energy evaluation                           
 ----------------------------------------------------------------------
```
以降に構造情報とエネルギーが以下のように出力されます。

```
Structure number  : 1
 File name         : ../xsf/structure01.xsf
 Number of atoms   : 2
 Number of species : 1
 
 Lattice vectors (Angstrom):
 
   a = (        1.40000200       0.00000000       0.00000000 )
   b = (        0.00000000       0.00000200       0.00000000 )
   c = (        0.00000000       0.00000000       0.00000200 )
 
 Cartesian atomic coordinates (input):
 
          x             y             z      
         (Ang)         (Ang)         (Ang)   
 --------------------------------------------
 C       0.000000     -0.000000     -0.000000
 C       1.400000     -0.000000     -0.000000
 
 Cohesive energy            :           -1.45972969 eV
 Total energy               :           -1.45972969 eV
```
ここで構造が``CRYSTAL``でなく``ATOMS``の場合、``Lattice vectors``は意味をなしません。
得られたエネルギーは解析的に計算して得られた値（``structure01.xsf``）-1.45972990 eVと良い一致を示していることが分かります。
最後に``04-predict/``で学習していない点（距離）についてエネルギーの予測を行います。
計算したい構造は``xsf/``に含まれています。予測のための入力ファイル（``predict.in``）は以下のようになります。

```predict.in
TYPES
1
C

NETWORKS
 C C.5t-5t.ann 2 5:tanh 5:tanh

FILES
301
./xsf/structure_001.xsf
./xsf/structure_002.xsf
./xsf/structure_003.xsf
./xsf/structure_004.xsf
./xsf/structure_005.xsf
./xsf/structure_006.xsf
./xsf/structure_007.xsf
./xsf/structure_008.xsf
./xsf/structure_009.xsf
./xsf/structure_010.xsf

...

./xsf/structure_301.xsf

```
``predict.x``の実行

```bash
predict.x predict.in > predict.out
```
得られた結果を図示すると以下のようになります。
![ljml.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/837941/aca25d9f-290e-f905-2052-773409c8abfa.png)

学習に用いた点が赤丸、予測した点を赤線で示しています。学習に用いた点では正確にエネルギーを再現している一方、それ以外の点では真の値と大きくずれていることが分かります。

## 水二量体
-ディレクトリ``Water_dimer``
ここでは水二量体のポテンシャルエネルギー曲線を学習してみましょう。この例ではMolnarらの構造を用いてrev-vdW-DF2汎関数を用いて計算したエネルギーを使用します。全エネルギーはQuantum-ESPRESSOを使用して計算しました。

``01-generate/``でトレーニングセットを生成します。

```generate.in
OUTPUT H2O.train

TYPES
2
H -13.5602905168717753
O -560.4742780438206798

SETUPS
H  H.fingerprint.stp
O  O.fingerprint.stp

FILES
14
../xsf/train4_newdimer_02.xsf
../xsf/train4_newdimer_03.xsf
../xsf/train4_newdimer_04.xsf
../xsf/train4_newdimer_05.xsf
../xsf/train4_newdimer_06.xsf
../xsf/train4_newdimer_07.xsf
../xsf/train4_newdimer_08.xsf
../xsf/train4_newdimer_09.xsf
../xsf/train4_newdimer_10.xsf
../xsf/train4_newdimer_11.xsf
../xsf/train4_newdimer_12.xsf
../xsf/train4_newdimer_13.xsf
../xsf/train4_newdimer_14.xsf
../xsf/train4_newdimer_15.xsf
```
この例では基底関数としてBehler-Parrinelloの対称関数を用います。フィンガープリントファイル（H.fingerprint.stp、O.fingerprint.stp）のパラメーターはMorawietz et al. [Proc. Natl. Acad. Sci. USA **113**, 8368 (2016)]のものを採用しましました。先の例と同様、``generate.x``を実行します。

```
generate.x generate.in > generate.out
```
これでトレーニングセットが``H2O.train``に書き出されます。
次に``02-train/``で学習を行います。

```train.in
TRAININGSET  H2O.train
TESTPERCENT  10
ITERATIONS   1000

TIMING

METHOD
bfgs

NETWORKS
 H  H.10t-10t.ann 2 10:tanh 10:tanh
 O  O.10t-10t.ann 2 10:tanh 10:tanh
```
``02-train/``で``01-generate/``にある``H2O.train``のシンボリックリンクを作ります。

```bash
ln -s ../01-generate/H2O`
```

次に``train.x``を実行します。

```bash
train.x train.in > train.out
```
``train.out``に

```train.out
 The optimization has converged. Training stopped.
 
 
 Training finished.
 
 ----------------------------------------------------------------------
                            Storing networks                           
 ----------------------------------------------------------------------
 
 Saving the H network to file : H.10t-10t.ann
 Saving the O network to file : O.10t-10t.ann
```
と表示されていれば最適化（トレーニング）は無事終了し、最終的なポテンシャルデータは``H.10t-10t.ann``と``O.10t-10t.ann``に出力されます。

今度は``03-test.in/``に移動し、ポテンシャルのテストを行います。先ずはポテンシャルへのシンボリックリンクを作成します。

```bash
ln -s ../02-train/H.10t-10t.ann
```

```bash
ln -s ../02-train/O.10t-10t.ann
```
以下の入力ファイルを準備します。

```predict.in
TYPES
2
H
O

NETWORKS
 H  H.10t-10t.ann 2 10:tanh 10:tanh
 O  O.10t-10t.ann 2 10:tanh 10:tanh

FORCES

FILES
17
../xsf/train4_newdimer_01.xsf
../xsf/train4_newdimer_02.xsf
../xsf/train4_newdimer_03.xsf
../xsf/train4_newdimer_04.xsf
../xsf/train4_newdimer_05.xsf
../xsf/train4_newdimer_06.xsf
../xsf/train4_newdimer_07.xsf
../xsf/train4_newdimer_08.xsf
../xsf/train4_newdimer_09.xsf
../xsf/train4_newdimer_10.xsf
../xsf/train4_newdimer_11.xsf
../xsf/train4_newdimer_12.xsf
../xsf/train4_newdimer_13.xsf
../xsf/train4_newdimer_14.xsf
../xsf/train4_newdimer_15.xsf
../xsf/train4_left.xsf
../xsf/train4_right.xsf
```

``predict.x``の実行

```bash
predict.in > predict.out
```
``predict.out``に出力された全エネルギーとDFT計算によって得られた全エネルギーを比較してみましょう。

最後に``04-predict``に移動し、``xsf/``以下にある構造のエネルギーを予測してみましょう。
ここでも先の作業と同様にポテンシャルへのシンボリックリンクを作成します。

```bash
ln -s ../02-train/H.10t-10t.ann
```

```bash
ln -s ../02-train/O.10t-10t.ann
```

以下のような入力ファイルを作成します。

```predict.in
TYPES
2
H
O

NETWORKS
 H  H.10t-10t.ann 2 10:tanh 10:tanh
 O  O.10t-10t.ann 2 10:tanh 10:tanh

FORCES

FILES
57
./xsf/train4_newdimer_01-02_000.xsf
./xsf/train4_newdimer_01-02_001.xsf
./xsf/train4_newdimer_01-02_002.xsf
./xsf/train4_newdimer_01-02_003.xsf
...
./xsf/train4_newdimer_14-15_001.xsf
./xsf/train4_newdimer_14-15_002.xsf
./xsf/train4_newdimer_14-15_003.xsf
./xsf/train4_newdimer_14-15_004.xsf
```
この例ではMolnarらの構造を元に水分子の間の重心間距離を内挿して新しい構造を作成しています。これまでと同様に

```bash
predict.x predict.in > predict.ou
```
を実行し全エネルギーを計算します。エネルギーのプロットのために``xsf/``以下に``dist_com``という重心間距離をÅで計算したデータがありますのでこれを用いて距離の関数としてエネルギーをプロットすると良いでしょう。

さてポテンシャルエネルギーの学習と予測はうまくいったのでしょうか？もしうまくいかないのであれば何が問題なのでしょうか？

## アルゴン二量体
-ディレクトリ``Ar_dimer``
この例ではrev-vdW-DF2を用いて計算したアルゴン二量体のエネルギー曲線を学習してみます。
先ず最初に``01-generate/``でトレーニングセットを生成します。

```generate.in
OUTPUT Ar2.train

TYPES
1
Ar -1287.322047978884

SETUPS
Ar Ar.fingerprint.stp

FILES
52
../xsf/ar2_3.00.xsf
../xsf/ar2_3.05.xsf
../xsf/ar2_3.10.xsf
...
../xsf/ar2_5.45.xsf
../xsf/ar2_5.50.xsf
../xsf/ar2_5.55.xsf
```
この例ではChebyshev多項式を使ったセットアップファイルを試してみます（レナードジョーンズ二量体と同様）。

```Ar.fingerprint.stp
DESCR
N. Artrith and A. Urban, Comput. Mater. Sci. 114 (2016) 135-150.
N. Artrith, A. Urban, and G. Ceder, Phys. Rev. B 96 (2017) 014112.
END DESCR

ATOM Ar

ENV 1
Ar

RMIN 1.0D0

BASIS type=Chebyshev
radial_Rc = 6.0 radial_N = 10
```
これまでと同様に``gerate.x``を実行します。

```bash
generate.x generate.in > generate.out
```
``02-train/``に移動してトレーニングセットのシンボリックリンクを作成します。

```bash
ln -s ../Ar2.train
```
以下の入力ファイルを使ってポテンシャルのトレーニングを行います。

```train.in
TRAININGSET Ar2.train
TESTPERCENT 20
ITERATIONS  1000

TIMING

METHOD
bfgs

NETWORKS
Ar Ar.5t-5t.ann 2 5:tanh 5:tanh
```
この計算では``TESTPERCENT``を10にすると最適化が1000回ではうまくいきませんでした。トレーニングセットの数を減らして``TESTPERCENT``を20にすることで数100回で最適化を行うことができました（ただこのあたりの動作がまだよく分かりません。うまくいっているのか、勝手にrestartしていて最適化にかかる回数が少なくなるのか...）。

``03-test/``に移動してポテンシャルへのシンボリックリンクの作成を行い、``predict.x``を実行します。

```bash
ln -s ../02-train/Ar.5t-5t.ann
```

```predict.in
TYPES
1
Ar

NETWORKS
Ar Ar.5t-5t.ann 2 5:tanh 5:tanh

FORCES

FILES
52
../xsf/ar2_3.00.xsf
../xsf/ar2_3.05.xsf
../xsf/ar2_3.10.xsf
...
../xsf/ar2_5.45.xsf
../xsf/ar2_5.50.xsf
../xsf/ar2_5.55.xsf
```

``predict.x``で計算されたエネルギーと力を、XSFファイルに記載されたDFTの値と比較してみましょう（単位に注意--referenceのXSFファイルの単位は面倒だったのでESPRESSOの出力のRy/Bohr単位のまま）。