AENETを使ってみる

はじめに

aenetを使ってポテンシャルエネルギー曲線の学習方法を学習してみます。
学習に用いたデータはこちらからダウンロードが可能です。

準備

aenetのコンパイル

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

mkdir -p ~/AENET/src

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を編集して環境に応じて適当なコンパイラとオプションを指定します。例えば

FC = ifort -c

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)]で紹介された結果を再現してみたいと思います。。

論文に記載されたパラメーターを使ってポテンシャルエネルギーを書いてみます。

対応する構造はダミー原子として炭素を使用して作成し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を実行します。

generate.x train.in > train.out

トレーニングセットはLJ.trainに出力されます。

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

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の実行

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の実行

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の実行

predict.x predict.in > predict.out

得られた結果を図示すると以下のようになります。

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

水二量体

-ディレクトリ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のシンボリックリンクを作ります。

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

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

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/に移動し、ポテンシャルのテストを行います。先ずはポテンシャルへのシンボリックリンクを作成します。

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

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の実行

predict.in > predict.out

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

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

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

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らの構造を元に水分子の間の重心間距離を内挿して新しい構造を作成しています。これまでと同様に

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を実行します。

generate.x generate.in > generate.out

02-train/に移動してトレーニングセットのシンボリックリンクを作成します。

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を実行します。

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単位のまま）。