GPU版LAMMPSによるテスト計算

# 概要
以前[WSL2上でCUDA環境を構築し，テストプログラムを実行する手続きについて解説した記事](https://qiita.com/jkoga/items/f7bd34e314ae23ebdfac)を書きました。本稿では，この時作った環境においてGPUで動作するプログラムを実際にビルドし，実行してみる試みを紹介したいと思います。対象のプログラムは，古典分子動力学シミュレーター[LAMMPS](https://www.lammps.org/)です。おおよそ以下のような内容を予定しています。
  
+ コンパイルする方法
    + GPUインターフェースライブラリーのビルド
    + LAMMPSロードモジュールのビルド 
+ 使い方
    + GPU上で走らせる方法
    + 複数GPUで実行する方法
+ テスト計算
    + 用いるポテンシャル EAMポテンシャル，TIP4Pポテンシャル
    + 環境 CPU : Intel core i7-12700 2.1 GHz, GPU : GeForce GTX 1650 super
+ むすび

用いるGPUはGPGPU専用製品と比較すると安価なグラフィックス描画用のものですが，このような環境で実アプリケーションにおいてどの程度の性能が見込めるのか確認していきたいと思います。

# コンパイルする方法

## GPUインターフェースライブラリー
LAMMPSインストールディレクトリーの下の `lib/gpu` にGPUのインターフェースライブラリーのソースコードが存在するので，これをビルドしていきます。

用いるGPUは倍精度の性能が単精度の1/32と低いため，GPUにおいてはできれば単精度で計算を実行したいです。そこで，「混合精度版」のインターフェースをビルドします。この版では，GPUの処理は単精度，CPUでの処理は倍精度で行います。最終的に重要なエネルギーや原子間力は倍精度で表すため，単精度を用いつつ計算精度の劣化を比較的小さくすることができます。

デフォルトが「混合精度版」なので，その実現のために特別なことを行う必要はありません。`Makefile.cuda_mps` において以下の行を編集することによって，このふるまいを制御することができます。
```
CUDA_PRECISION = -D_SINGLE_DOUBLE
```
`_SINGLE_DOUBLE` のかわりに `_SINGLE_SINGLE` とすればすべて単精度， `_DOUBLE_DOUBLE` とすればすべて倍精度で計算を行うようになります。

以前作った環境において以下のコマンドを実行し，[NVIDIA HPC SDK](https://developer.nvidia.com/hpc-sdk)をMPI込みで有効にします。
```
module load nvhpc-hpcx-cuda12
Loading nvhpc-hpcx-cuda12/23.9
  Loading requirement: hpcx
```
用いるコンパイル環境はNVIDIA HPC SDKなので，それに対応した `Makefile.cuda_mps` という`Makefile`を用います。すなわち，以下の要領でビルドします。
```
cd LAMMPS_INSTALL_DIR/lib/gpu
make -f Makefile.cuda_mps
nvcc -I/usr/local/cuda/include -DUNIX -O3 --use_fast_math -DLAMMPS_SMALLBIG -Xcompiler -fPIC -gencode arch=compute_50,code=[sm_50,compute_50] -gencode arch=compute_52,code=[sm_52,compute_52] -gencode arch=compute_60,code=[sm_60,compute_60] -gencode arch=compute_61,code=[sm_61,compute_61] -gencode arch=compute_70,code=[sm_70,compute_70] -gencode arch=compute_75,code=[sm_75,compute_75] -gencode arch=compute_80,code=[sm_80,compute_80] -gencode arch=compute_86,code=[sm_86,compute_86] -D_SINGLE_DOUBLE --fatbin -DNV_KERNEL -o .//amoeba.cubin .//lal_amoeba.cu
...
...
```
このようなコマンドによって`libgpu.a` というライブラリーが生成されれれば成功です。

## LAMMPS本体のビルド

LAMMPS本体をコンパイルし，`libgpu.a` にリンクしてロードモジュールを作成します。そのためパッケージ `gpu` を以下の要領で有効にします。その他，テスト計算で用いる予定のポテンシャルのパッケージを有効にしています。作業はLAMMPSインストールディレクトリーの下の`src`ディレクトリーにおいて行います。
```
cd ../../src
make yes-gpu
make yes-manybody
make yes-kspace
make yes-rigid
make mpi
Gathering installed package information (may take a little while)
...
...
```

ロードモジュール`lmp_mpi` が生成されれば無事コンパイルできたことになります。

# 使い方
コンパイルしたロードモジュールをGPU上で実行するには二つの方法があります。
+ ポテンシャル指定の最後に `/gpu` をつける。例： `lj/cut/gpu`
+ `-sf gpu` をつけて実行する。例： `mpirun -np 2 lmp_mpi -in in -sf gpu`；ここで `-in in`はインプットファイル指定

複数のGPUが搭載されているようなシステムの場合，以下の要領で複数のGPU上で動作させることも可能です。
```
mpirun -np 2 lmp_mpi -in in -sf gpu -pk gpu 2
```
`-pk gpu 2`によって2基のGPUを用いる設定になります。多くの場合MPIプロセス数とGPU数を一致させた方が良好なパフォーマンスを得ることができると思いますが，複数のGPUを搭載したマシンを使っている場合は試行錯誤して最良の組み合わせを見出すのがよいのではないかと思います。

# テスト計算
以下のような問題を用いてテスト計算を実行してみました。
+ EAMポテンシャルを用いたbcc-Feスーパーセル ($100 \times 100\times 100$) , 200万原子の20ステップのMD計算
+ TIP4Pポテンシャルを用いた2484分子の1000ステップのMD計算

計算負荷の観点からは，二つの例題は長距離相互作用が含まれるかどうかという点が最も大きな違いです。すなわち，EAMポテンシャルは短距離相互作用のみなのに対し，TIP4Pポテンシャルは長距離相互作用であるクーロン相互作用を含みます。クーロン相互作用はPPPM法によって評価します。なお，[LAMMPSのマニュアル](https://docs.lammps.org/kspace_style.html)によると汎用のPPPM法は一部の演算をGPUにおいて計算することが可能なようなのですが，TIP4Pで用いる専用PPPM法 (`pppm/tip4p`) はGPU非対応なのでCPUのみで動作します。

動作環境は以下の通りです。
+ CPU : [Intel core i7-12700 2.1 GHz](https://www.intel.com/content/www/us/en/products/sku/134591/intel-core-i712700-processor-25m-cache-up-to-4-90-ghz/specifications.html)
+ GPU : [GeForce GTX 1650 super](https://www.nvidia.com/ja-jp/geforce/graphics-cards/16-series/)

用いたCPUには8つのperformance coreが搭載されています。そこで，8並列までは性能向上が期待できると考えMPIプロセス数は非並列と8並列にしました。GPUは1基しか搭載されてないので，MPI非並列実行のみです。

計算に要した時間は以下の通り（単位は秒）。

|     | EAM bccFe 200万原子 | TIP4P water 2484分子 |
|-----|--------------------|---------------------|
| GPU | 1.9                | 7.0                 |
| CPU 非並列 | 23.2         | 37.7 |
| CPU 8並列 | 6.4 |9.9 |

表の数値から明らかなように，GPU上でCPUよりも高速に動作することが確認できました。
+ EAMポテンシャルの場合は非並列と比較して20倍以上，8並列と比較しても3倍以上の高速化が達成できています
+ TIP4Pポテンシャルの場合差はそこまで大きくなく，非並列の場合は5倍以上，8並列の場合は1.4倍となりました

TIP4Pポテンシャルの場合相対的にGPUが低速な理由は，PPPM法によるクーロン相互作用の計算がGPUに対応していないことです。GPU化されたPPPM法を用いることができれば改善する可能性があります。ただし[マニュアルによると](https://docs.lammps.org/Speed_gpu.html)GPU化がなされているのは一部で，アルゴリズム上必須の高速フーリエ変換はCPUで行われるとのことなので，どこまで効果があるかは試してみないと分かりません。

LAMMPSは，ログファイルの最後の方でどの演算にどのくらいの時間がかかったかを記したレポートを出力してくれます。TIP4Pの場合にその内容を紹介したいと思います。

GPU：
```
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0.83107    | 0.83107    | 0.83107    |   0.0 | 11.90
Bond    | 0.12234    | 0.12234    | 0.12234    |   0.0 |  1.75
Kspace  | 3.7883     | 3.7883     | 3.7883     |   0.0 | 54.25
Neigh   | 0.00084586 | 0.00084586 | 0.00084586 |   0.0 |  0.01
Comm    | 0.075925   | 0.075925   | 0.075925   |   0.0 |  1.09
Output  | 0.0044897  | 0.0044897  | 0.0044897  |   0.0 |  0.06
Modify  | 2.1266     | 2.1266     | 2.1266     |   0.0 | 30.46
Other   |            | 0.0328     |            |       |  0.47
```

CPU 非並列：
```
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 28.519     | 28.519     | 28.519     |   0.0 | 75.66
Bond    | 0.11389    | 0.11389    | 0.11389    |   0.0 |  0.30
Kspace  | 3.8728     | 3.8728     | 3.8728     |   0.0 | 10.27
Neigh   | 2.891      | 2.891      | 2.891      |   0.0 |  7.67
Comm    | 0.079032   | 0.079032   | 0.079032   |   0.0 |  0.21
Output  | 0.004175   | 0.004175   | 0.004175   |   0.0 |  0.01
Modify  | 2.1671     | 2.1671     | 2.1671     |   0.0 |  5.75
Other   |            | 0.04474    |            |       |  0.12
```

CPU 8並列：
```
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 4.801      | 5.0749     | 5.4564     |   8.2 | 53.23
Bond    | 0.018259   | 0.019488   | 0.02023    |   0.4 |  0.20
Kspace  | 1.9346     | 2.3158     | 2.5893     |  12.1 | 24.29
Neigh   | 0.56025    | 0.56047    | 0.56073    |   0.0 |  5.88
Comm    | 0.15301    | 0.15364    | 0.15407    |   0.1 |  1.61
Output  | 0.0017384  | 0.0019381  | 0.0026595  |   0.6 |  0.02
Modify  | 1.2191     | 1.2713     | 1.302      |   2.4 | 13.33
Other   |            | 0.1373     |            |       |  1.44
```

この結果について何点かコメントしたいと思います。
+ GPUでは`Pair` および `Neigh` の計算時間が際立って短いです。その他の部分はCPU非並列の場合とほぼ同等で，GPU化がなされていないと考えられます
+ CPU非並列と8並列を比較すると，`Pair`, `Neigh` などの演算は5～6倍程度高速化されています。`Kspace` も速くなってはいますが，上であげた二つほどではなく2倍程度です。これは高効率の並列化が比較的難しい高速フーリエ変換を含む処理だからではないかと考えられます

# むすび
本稿では，以前構築したGPU環境を用いて古典分子動力学シミュレーターLAMMPSによるテスト計算を行ってみました。コンパイルや実行の仕方を説明し，EAMポテンシャルを用いたbccFeとTIP4Pポテンシャルを用いた水のシミュレーションの計算時間を計測してみました。その結果，グラフィックス処理用でミドルクラスの比較的安価なGPUでも高速な計算を行えることが確認できました。特に，長距離相互作用を伴わない場合に大きな効果を発揮しました。今後機会があれば[GPGPU専用製品](https://ja.wikipedia.org/wiki/NVIDIA_Tesla)を用いた検証も行っていきたいと思います。