Bug: some inter atomic pair forces are not correctly calculated when number of pairs is large

[passive-radio]

Fact

1. 一部の atomic pair でCURPの計算ができていない or できているが出力時に問題がある

Tested environment

• CURP v1.3.1
• CURP v2

[passive-radio]

Cited from lab's internal mattermost post.

計算環境 windows11+vscodeでBassへssh接続して実行 curp v1.3.1, openmpi 4.1.5 計算機を用いた4コア並列計算 使ったスクリプト インプット群, jobサブミット用スクリプト, fluxのジョブスクリプト, conductivityのジョブスクリプト

とりうるすべての二次構造グループのペアを作り、計算を実行した サブミット用スクリプトがあるところでbash submit_inter_secondary.sh を実行すると、startからendまでの番号のトラジェクトリそれぞれについてflux、conductivityを順番に計算してくれる バグの内容 上のスクリプトをbass上で計算すると、一部のペアでconductivityがnanとなる(outputをこのスレッドのリプライにぶら下げておきます) nanとなるペアはいずれも一次構造的に連続していない(bonded interactionを含まない)なグループのペア nanとなるペアのトラジェクトリによる差異はない、つまり、すべてのトラジェクトリについて特定のペアで計算結果がnanとなる これまでに分かっていること 原因はfluxの計算にあるっぽい flux.ncをのぞくとconductivityの結果がnanとなっているペアのfluxの値がすべてのtime stepで0になっていた ペア数を極端に減らすとnanだったペアの計算が正常に実行されているように見える(アウトプットをぶら下げておきます) curpのversionを1.3.1->2.0にしても同様の現象がみられる 並列計算を1スレッドにしても結果は変わらない

[passive-radio]

問題に関連する処理の流れ。以下のどこか、または複数でバグが存在する可能性が非常に高い

1. https://github.com/yamatolab/current-calculations-for-proteins/blob/b27ab52095143429e264b0b8224cbcdfe1144410/curp/table/interact_table.py#L6 がどの i, j atom で２体力を計算するかを決める interaction table list を作成する
2. 作成された i,j interactions を定義する interaction_table list が https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/compute.py#L292 で Python で実装された Calculator クラスに渡される。
3. flux, current calculator はいずれも CalculatorBase 基底クラスを継承していて、いずれの calculator クラスも CalculatorBase.prepare メソッドをそのまま使っている。 https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/current/base.py#L33 ここで、受け取った interaction table list を TwoBodyCalculator の実体に渡している。
4. interact table liist を受け取っている TwoBodyCalculator は https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/twobody/__init__.py#L7 で force field が amber系だったときに https://github.com/yamatolab/current-calculations-for-proteins/blob/v1.3.1/curp/twobody/amberbase.py で定義された TwoBodyForceBase クラスのこと。
5. 受け取った interact table list を使って https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/twobody/amberbase.py#L60 で TwoBodyForceBase.cal_force メソッドが、coulomb, vdw を計算するメソッドに対して interact table list の各子要素 interact table を渡して、２体力を計算させている
6. cal_coulomb, cal_vdw メソッドの実際の処理は _cal_nonbonde メソッド https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/twobody/amberbase.py#L201 になる。
7. _cal_nonbonded メソッドは、https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/twobody/amberbase.py#L213 の通り fortran で定義された２体力計算 module を使って mod.calculate(table) しており、受け取った interaction table を元に２体力を計算している
8. interaction table ごとに計算された２体力などの値が https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/twobody/amberbase.py#L235 で "energy", "forces", "tbforces", "displacement" がキーの dict の値として格納されて _cal_nonbonded メソッド (つまり、それを実行している cal_coulomb, cal_vdw メソッド）などで return される。
9. この cal_coulomb, cal_vdw メソッドは https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/current/flux.py#L231 や https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/current/flux.py#L253 などの FluxCalculator 内でジェネレーターとして呼ばれて、https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/current/flux.py#L257 などで interact table list の子要素 (小分けにされた interact table) ごとにジェネレーターが実行されて "tbforces" キーの値を取得することで２体力を取得して、flux 計算用のクラスに与えている。
10. そして https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/current/flux.py#L298 などで受け取った２体力をもとに実際に flux を計算して https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/current/flux.py#L328 で flux の計算結果を返している。

[passive-radio]

以上が熱流を計算するペアを指定した interaction 設定ファイルを読み込んでから、実際に２体力を計算し、熱流を計算して python 側が受け取るまでの流れ。

[passive-radio]

current-calculations-for-proteins/curp/table/interact_table.py

Line 6 in b27ab52

class InteractionTableGenerator: がどの i, j atom で２体力を計算するかを決める interaction table list を作成する

このあたりがまず怪しそう。

このクラスでやっていることは、 計算したい残基ペアや２次構造のペアなどの interactions を定義した設定ファイルから fortran の２体力計算モジュールに渡す interction table を作成して渡している。 だが、一度に大量の interaction を計算させるとメモリに乗り切らずエラーになってしまう。なので、 fortran 計算モジュールで計算させる atom i,j のリスト (interaction table) を小分けにして、渡す機能が実装されている。その肝心の interaction がある一定数を超えると小分けにする実装が以下にされている。

https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/table/interact_table.py#L418

ここでは、一度に計算させる interaction の数を https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/table/interact_table.py#L196 で 10 と定義した _memory_limit を使って、 https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/table/interact_table.py#L427

のように _memory_limit * 10**6 * 1.2 と定義している。なので、例えば計算させたい atom i,j のinteractions の数が ２次構造間の熱流のように非常に大きい場合、10**7 pairs を超える場合がある。 そうすると、計算させた interaction 全てを一度に計算させるのではなく interaction table を小分けにする処理が

https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/table/interact_table.py#L441

https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/table/interact_table.py#L447

で実行される。

なぜ怪しいのかというと、一度に渡す interractions 数にキャップを付ける機能を司どると言っても良い、https://github.com/yamatolab/current-calculations-for-proteins/blob/76804fd3d612e113815289cabc5a8821c2285457/curp/table/interact_table.py#L196 の _memory_limit を既存の 10 から 1000 にすると、正常に計算結果が返ってきたという報告があるため。

つまり、interaction table の小分けが行われると fortran の計算 or 計算結果を python が受け取る or 出力する箇所でバグが発生している可能性が非常に高い。

なので、interaction table が小分けになった場合とならない場合で処理にどのような差異が出るかを追っていけばバグの原因がわかりそう.