AB Initio分子動力学シミュレーションからの溶融物と流体の分析

1. 分子動力学の分析 注: このパッケージは、GitHub のウェブサイト (https://github.com/rcaracas/UMD_package) から、オープンアクセス パッケージとして ERC IMPACT プロジェクトの専用ページ (http://moonimpact.eu/umd-package/) を介して入手できます。 パッケージから 1 つ以上の専用 Python スクリプトを使用して、特定の物理プロパティのセットを抽出します。コマンド ラインですべてのスクリプトを実行します。これらはすべて一連のフラグを採用しており、スクリプト間で可能な限り一貫性があります。フラグ、その意味、およびデフォルト値はすべて 表 1 にまとめられています。 旗 意味 それを使用してスクリプト 既定値 -h 短いヘルプ すべての -f UMD ファイル名 すべての -i 廃棄する熱化ステップ すべての 0 -i 原子間結合を含む入力ファイル 種 分化 債券.インプット -s 周波数のサンプリング msd, 種分化 1 (すべてのステップが考慮されます) -a 原子または陰イオンの一覧 種 分化 -c カチの一覧 種 分化 -l ボンド長さ 種 分化 2 -t 温度 振動,レロジー -v 平均二乗変位解析のための軌道のサンプリングウィンドウの幅の分離 msd 20 -z 平均二乗変位解析のための軌道のサンプリングウィンドウの開始の分離 msd 20 表 1: UMD パッケージで使用される最も一般的なフラグとその最も一般的な重要性。 まず、VASP8 や QBox9 などの第一原理コードで実行される MD シミュレーションの出力を UMD ファイルに変換します。 MD シミュレーションが VASP で実行された場合は、コマンド ラインで次のように入力します。VaspParser.py -f -i ここで–f フラグは VASP OUTCAR ファイルの名前を定義し、–i の熱化の長さを定義します。注: –i で定義された最初のステップでは、シミュレーションの最初のステップを破棄できます。典型的な分子動力学の実行では、計算の最初の部分は、温度の温度化、すなわち、すべての原子が温度のガウス状分布を表すためにシステムを取り、平衡値の周りに温度、圧力、エネルギーなどの変動を示すためにシステム全体のためにかかる時間を表します。流体の統計的特性を解析する際には、シミュレーションのこの熱化部分を考慮しないでください。 を変換します。umd ファイルを に変換します。VMD4 や Vesta5 など、さまざまなパッケージで視覚化を容易にする xyz ファイル。コマンド ラインで次のように入力します。umd2xyz.py -f -i -s ここで –f は の名前を定義します。umd ファイル -i は破棄する熱化期間を定義し、–s に格納されている軌道のサンプリングの頻度を定義します。umd ファイル。既定値は –i 0 –s 1、つまり、破棄されずにシミュレーションのすべてのステップを考慮します。 umd2poscar.py スクリプトを使用して 、umd ファイルを VASP タイプの POSCAR ファイルに逆にします。シミュレーションのスナップショットは、事前に定義された頻度で選択できます。コマンド ラインで次のように入力します。umd2poscar.py -f -i -l -s ここで –l は POSCAR ファイルに変換される最後のステップを表します。デフォルト値は -i 0 -l 1000000 -s 1 です。–l のこの値は、一般的な全体の軌道をカバーするのに十分な大きさです。 2. 構造解析を実行する gofrs_umd.pyスクリプトを実行して、アトミックタイプAとBのすべてのペアのペア分布関数 (PDF) g ᴀʙ(r) を計算します (図 3)。出力は、タブ区切りの 1 つの ASCII ファイルに、拡張 gofrs.dat 書き込まれます。コマンド ラインで次のように入力します。gofrs_umd.py -f -s -d -i 注: デフォルトはSampling_Frequency(軌道をサンプリングする周波数)=1ステップです。分離間隔間隔 (g(r)) = 0.01 Å;InitialStep (破棄される軌道の先頭のステップ数) = 0。放射状PDF、g ᴀʙ(r)は、A型の原子を中心とした半径rと厚さdrの球状シェル内の距離d_ᴀʙ B型の平均数です(図3)。ρでは、原子密度、NAおよびNBタイプAおよびBの原子数、およびδ(r-r ᴀʙ)は、原子AとBがrとr+drの間の距離にある場合に1に等しいデルタ関数を示します。最初の最大値の g ᴀʙ(r) の棄権は、A 型と B 型の原子の間で最も高い確率の結合長さを示します。最初の最小値は、最初の座標球の範囲を区切ります。したがって、最初の最小値までのPDF上の積分は、平均調整数を与えます。A型とB型のすべての組み合わせに対するg ᴀʙ(r)のフーリエ変換の合計は、回折計で実験的に得られたように、流体の回折パターンを生み出す。しかし実際には、高次の配位球がg ᴀʙ(r)から欠落している場合、回折パターン全体を得ることができない場合があります。 図3:ペア分布関数の決定(a) 1つの種の各原子(例えば赤)について、配位種のすべての原子(例えば灰色および/または赤)は距離の関数として数えられる。(b) 各スナップショットの結果として得られる距離分布グラフは、この段階ではデルタ関数の集まりに過ぎず、すべての原子とすべてのスナップショットを平均化し、(c)連続するペア分布関数を生成する理想的なガス分布によって重み付けされる。 最初の最小値は 、最初の配位球の半径で、後で種分分析で使用します。 この図の大きなバージョンを表示するには、ここをクリックしてください。 最初の配位球の半径として平均原子間結合距離を抽出します。このためには、 g ᴀʙ(r)関数の最初の最大値の位置を特定します: gofrs.datファイルをスプレッドシートアプリケーションでプロットし、原子の各ペアの最大値と最小値を検索します。 スプレッドシート ソフトウェアを使用して、最初の調整球の半径を PDF の最初の最小値である g ᴀʙ(r)として識別します。これは流体の構造解析全体の基礎となります。PDF は流体中の原子の平均結合状態を生成します。 最初のミニマ、すなわち、アブシッサの距離を抽出し、例えば bonds.inputと呼ばれる別のファイルに書き込みます。または、umD パッケージの analyze_gofr スクリプトのいずれかを実行して 、g ᴀʙ(r) 関数の最大値と最小値を識別します。コマンド ラインで次のように入力します。analyze_gofr_semi_automatic.py プログラムで開いたグラフに表示される g ᴀʙ(r) 関数の最大値と最小値の位置をクリックします。スクリプトは、現在のフォルダを自動的にスキャンし、すべての gofrs.datファイルを識別し、それぞれのファイルの分析を実行します。スクリプトが教育的な初期推測を必要とするたびに、ウィンドウの最大値と最小値をもう一度クリックします。 開いて、原子間の結合距離を含む bonds.input と呼ばれる自動的に生成されたファイルを見ます。 3. 種分分析を実行する グラフ理論内の接続性の概念を使用して、原子間の結合のトポロジを計算します:原子はノードであり、原子間結合はパスです。 speciation_umd.py スクリプトには bonds.input ファイルで定義された原子間結合距離が必要です。注: 接続マトリックスは、各時間ステップで構築されます: 対応する最初の配位球の半径よりも小さい距離にある 2 つの原子は、結合されていると見なされます。つまり接続されています。さまざまな原子ネットワークは、この幾何学的基準によって接続が定義されたグラフ内のノードとして原子を扱うことによって構築されます。これらのネットワークは原子種であり、そのアンサンブルは、その特定の流体における原子種分化を定義する(図4)。 図4:原子クラスタの識別配位多面体は、原子間距離を使用して定義されます。指定した半径よりも小さい距離にあるすべての原子は、結合されていると見なされます。ここでの閾値は、図1で定義された第1の配位球(明るい赤い円)に相当する。重合および従って化学種は、結合原子のネットワークから得られる。無限のポリマーを形成する他の原子から分離された中央のRed1Grey2クラスターに注意してください。この図の大きなバージョンを表示するには、ここをクリックしてください。 スペシエーション スクリプトを実行して接続行列を取得し、配位多面体または重合を取得します。コマンド ラインで次のように入力します。speciation_umd.py -f -s -i -l -c -a -m -r ここで-i フラグは、前のステップで生成された、例えば、原子間結合距離を持つファイルを示します。または、-l フラグで定義されたすべての結合に対して 1 つの長さのスクリプトを実行します。注: -c フラグは、中心の原子を指定し、-a フラグをリガンドに指定します。中央原子とリガンドの両方が異なるタイプであることができます。この場合、カンマで区切る必要があります。m フラグは、分析で考慮されるために種が生きなければならない最小時間を示します。デフォルトでは、この最小時間はゼロで、最終的な分析でカウントされるすべてのオカレンスです。 第 1 レベルで接続グラフをサンプリングして、コーディネーション多面体を識別するフラグ –r 0 を指定してspeciation_umd.pyスクリプトを実行します。例えば、カ チオン として示される中心原子は、1つ以上 の陰イオン で囲まれていてもよい(図4)。スペシレーションスクリプトは、配位多面体の一つ一つを識別します。すべての配位多面体の加重平均は、調整数を、PDFの統合から得られたものと同一に与える。コマンド ラインで次のように入力します。speciation_umd.py -f -i -c -a -r 0注: 流体の平均配位数は小数です。この分数は、コーディネートの平均特性から生まれます。種分化に基づく定義は、流体の構造をより直感的で有益な表現をもたらす。 speciation_umd.py スクリプトを 実行するには、フラグ –r 1 を指定して、接続性グラフをすべての深度レベルでサンプリングして重合を取得します。原子グラフを通るネットワークは、原子が他の結合(例えば、交互のカチオンおよび陰イオンの配列)にさらに離れて結合されるので、ある深さを有する(図4)。 2 つのファイルを開きます。ポパル.dat と.統計.dat 連続して;これらは、スペシテーション スクリプトの出力を構成します。各クラスターは、その化学式、それが形成された時間、それが死んだ時間、その寿命、このクラスターを形成する原子のリストを持つ行列を指定して、1行に書かれています。.popul.dat ファイルに示されているように、シミュレーションで見つかったすべての化学種の各原子クラスターの寿命をプロットします (図 5)。 で見られる、各種の豊富な人口分析をプロットします。stat.dat ファイル。この分析は、絶対と相対の両方で、ケース -r 0 の座標多面体の実際の統計量に対応します。重合の場合、-r1では、原子の相対数に対する正規化を適用する必要がある場合があるため、慎重に扱う必要があります。豊富さは、生涯にわたって積分に対応します。ザ。stat.dat ファイルには、各クラスターのサイズ、すなわち、そのクラスターを形成する原子の数もリストされます。 4. 拡散係数の計算 自己拡散性を得るために時間の関数として原子の平均二乗変位(MSD)を抽出する。MSD の標準式は次のとおりです。ここで、プリファクターは再正規化です。MSD ツールを使用すると、流体の動的な側面を分析するさまざまな方法があります。注: T はシミュレーションの合計時間であり、N αはタイプ αの原子数です。初期時間 t0 は任意であり、シミュレーションの前半に及びます。ニニトは初期回数です。τ は、MSD が計算される時間間隔の幅です。最大値はシミュレーションの時間の長さの半分です。一般的な MSD 実装では、すべてのウィンドウは前のウィンドウの最後から開始されます。しかし、まばらなサンプリングは、MSDの傾きを変更することなく、MSDの計算を高速化することができます。このため、i 番目のウィンドウは時刻 t0(i) から始まりますが、(i+1)番目のウィンドウは時刻 t0(i) + τ + v で始まり、v の値はユーザー定義です。同様に、ウィンドウの幅は、ユーザーによって定義された個別のステップで増加します。z (“水平ステップ”) と v (“垂直ステップ”) の値は正またはゼロです。両方のデフォルトは 20 です。 一連の msd_umd スクリプトを使用して MSD を計算します。出力は.msd.dat ファイルでは、各原子タイプ、原子、またはクラスタの MSD が時間の関数として 1 つの列に出力されます。 各アトミックタイプの平均 MSD を計算します。MSD は各原子に対して計算され、その後、アトミックタイプごとに平均化されます。出力ファイルには、アトミック・タイプごとに 1 つの列が含まれます。コマンド ラインで次のように入力します。msd_umd.py -f -z -v -b 各原子の MSD を計算します。MSD は各原子に対して計算され、その後、アトミックタイプごとに平均化されます。出力ファイルには、シミュレーション内のアトムごとに 1 つの列が含まれ、次にアトミックタイプごとに 1 つの列が含まれます。この機能により、液体や気体、または2つの液体など、2つの異なる環境で拡散する原子を識別できます。コマンド ラインで次のように入力します。msd_all_umd.py -f -z -v -b 化学種の MSD を計算します。スペシテーション スクリプトで識別され、 で印刷されるクラスターの母集団を使用します。ポプル.dat ファイル。MSD は、個々のクラスタごとに計算されます。出力ファイルには、クラスターごとに 1 つの列が含まれます。大規模なポリマーの検討を避けるために、クラスターのサイズに制限を設けます。デフォルトは 20 個の原子です。コマンド ラインで次のように入力します。msd_cluster_umd.py -f -p -s -b -c 注: デフォルト値は–b 100 –s 1 –c 20 です。 スプレッドシートベースのソフトウェアを使用して MSD をプロットします(図 6)。MSD 対時間のログ ログ表現で、勾配の変化を特定します。弾道 体制を 表す最初の部分、通常は短い、すなわち衝突後の原子の速度の保全を分離する。2番目の長い部分は、 拡散 体制、すなわち衝突後の原子の速度の散乱を表します。 MSD の傾きから拡散係数を次のように計算します。ここで、Z は自由度の数(平面内の拡散の場合は Z = 2、空間での拡散の場合は Z = 3)、t は時間ステップです。 5. 時間相関関数 一般式を使用して、システムの慣性の尺度として時間相関関数を計算します。A は、原子の位置、原子速度、応力、偏光などのさまざまな時間依存変数であり、それぞれの降伏は、緑久保関係12,13を介して、異なる物理的特性を介して、時にはさらに変換された後に起きます。 原子速度を解析し、液体の振動スペクトルと原子自己拡散係数の代替式を得る。 vibr_spectrum_umd.pyスクリプトを実行して、各アトミックタイプの原子速度自動相関(VAC)関数を計算し、高速フーリエ変換を実行します。コマンド ラインで次のように入力します。vibr_spectrum_umd.py -f -t ここで –t は、ユーザーが定義する必要がある温度です。このスクリプトは、2 つのファイルを出力します。vels.scf.dat 各アトミックタイプの VAC 関数を持つファイルと.vibr.dat 各原子種で分解された振動スペクトルと合計値を持つファイルです。 を開いて、 vels.scf.datを読みます。スプレッドシートのようなソフトウェアを使用して vels.scf.dat ファイルからVAC関数をプロットします。 フーリエVACの本当の部分を維持します。これは、周波数の関数として、振動スペクトルを生成するものです。ここで m は原子の質量です。 スプレッドシートのようなソフトウェアを使用して 、vibr.dat ファイルから振動スペクトルをプロットします(図7)。流体のディフューシャブ特性と有限周波数のスペクトルのさまざまなピークに対応する 、ω=0の有限値を特定します。振動スペクトルへの各原子型の参加を特定します。注: 原子タイプの分解は、異なる原子が拡散係数に対応する 異なるω=0寄与を有することを示しています。スペクトルの一般的な形状は、対応するソリッドよりも少ないフィーチャではるかに滑らかです。 シェルでは、各原子種の拡散係数をもたらす振動スペクトルの積分を読み取ります。注:熱力学的特性は、振動スペクトルからの統合によって得ることができますが、結果は2つの近似のために注意して使用する必要があります:統合は、必ずしも高温で保持されない準高調波近似の範囲内で有効です。拡散に対応するスペクトルの気体状部分を廃棄する必要があります。統合は、スペクトルの格子状の部分でのみ行われるべきです。ただし、この分離には通常、さらにいくつかの処理後の手順と計算が必要です14、これは現在の UMD パッケージではカバーされません。 viscosity_umd.pyスクリプトを実行して、成分応力テンソルの自己相関を分析して、溶融の粘度を推定します。コマンド ラインで次のように入力します。viscosity_umd.py -f -i -s -o -l 注: この機能は探索的なものであり、結果は慎重に行う必要があります。まず、シミュレーションの長さに対する粘度の収束を十分に確認します。 ストレステンソル15 の自己相関から流体の粘度を次のように導出します。ここで、VとTはそれぞれ体積と温度で、κBはボルツマン定数であり、σ ijはデカルト座標で表される応力テンソルのij対角外成分である。 より適切な適合度を使用して、粘度15,16のより強い推定値を得て、有限サイズとシミュレーションの有限期間から生じる可能性のある応力-テンソル自動相関関数のノイズを回避します。応力テンソルの自動相関関数の場合は、良好な結果をもたらす次の機能 form15,16 を使用します。ここで、A、B、τ1、τ2、およびωは適合パラメータです。 統合後、粘度の式は次のようになります。 6.シミュレーションに由来する熱力学パラメータ。 averages.py 実行して、圧力、温度、密度、および内部エネルギーの平均値と広がり(標準偏差)を umd ファイルから抽出します。コマンド ラインで次のように入力します。averages.py -f -s の値は 、デフォルトとして –s 0 を使用します。 ブロッキングメソッドを使用して、平均の統計的誤差を計算します。注:この方法にはさまざまな種類があります。アレンとティルデスリー2の仕事に続いて、時間ブロックのシーケンスにわたって平均を平均し、ますます長い、および算術平均17に関して標準偏差を推定することが一般的です。サンプリングが相関していない場合、収束は、多くて長いブロック サイズの限界に達する可能性があります。収束の実際のしきい値は通常手動で選択する必要があります。 halving method18を使用する:最初のデータサンプルから始まり、各ステップ のκで、前のステップ のκ-1から2つの対応する連続したサンプルごとに平均化することによってサンプルの数を半分にします。 fullaverages.py スクリプトを実行して、平均の誤差を含む完全な統計分析を実行します。コマンド ラインで次のように入力します。fullaverages.py -s -u 注: スクリプトは、現在のディレクトリ内のすべての .umd.dat ファイルを検索し、それらのすべての分析を実行する時点に自動化されています。デフォルトは –s 0 –u 0 です。u 0 の場合、出力は最小限であり、-u 1 の出力はフルであり、いくつかの代替単位が出力されます。このスクリプトは、平均値の誤差を推定するために収束をチェックするためのグラフィカルイメージを作成するため、グラフィカルなサポートが必要です。