The Analysis of Dissociation Process of Carbonic Acid in Aqueous by Using Molecular Dynamics with Machine Learning Potential

Abstract

For example， dissolved carbon dioxide and carbonic acid in aqueous solution are important compounds for life activities， serving as a buffer system for pH regulation in the body. Theoretical calculations for carbon dioxide and carbonic acid， and the transition states， which are the states mediated by water molecules have been vigorously performed， however， there are limitations in terms of computational speed and resources for exploring various conformations and reaction mechanisms in the water molecule network through ab initio calculations. The machine learning potential (MLP) was created by using the DeepMD framework， which is one of the MLPs that have been reported to be applied in various systems. We analyzed the conformation of the dissociation process of carbonic acid by using the Nudged Elastic Band method in the molecular dynamics (MD) calculations. The energy difference between MLP-MD and ab-initio MD (AIMD) at the beginning and the end of the reaction was 0.109 kcal/mol/atoms. Looking at the trend of the barrier value with respect to the number of water molecules in the water molecular network， the barrier value to decrease with the increase in the number of water molecules has similar tendency of the previous quantum chemical and AIMD calculations.

Translated Abstract

For example， dissolved carbon dioxide and carbonic acid in aqueous solution are important compounds for life activities， serving as a buffer system for pH regulation in the body. Theoretical calculations for carbon dioxide and carbonic acid， and the transition states， which are the states mediated by water molecules have been vigorously performed， however， there are limitations in terms of computational speed and resources for exploring various conformations and reaction mechanisms in the water molecule network through ab initio calculations. The machine learning potential (MLP) was created by using the DeepMD framework， which is one of the MLPs that have been reported to be applied in various systems. We analyzed the conformation of the dissociation process of carbonic acid by using the Nudged Elastic Band method in the molecular dynamics (MD) calculations. The energy difference between MLP-MD and ab-initio MD (AIMD) at the beginning and the end of the reaction was 0.109 kcal/mol/atoms. Looking at the trend of the barrier value with respect to the number of water molecules in the water molecular network， the barrier value to decrease with the increase in the number of water molecules has similar tendency of the previous quantum chemical and AIMD calculations.

1 Introduction

生体内のpH調整の緩衝系として，水溶液中の溶存二酸化炭素(CO₂)と炭酸(H₂CO₃)は生命活動にとって重要な化合物である．生体分子の機能解析にとって水ネットワーク内中のCO₂，重炭酸イオン，H₂CO₃の定量的な解析は，新機能を見出すための基盤的な知見を与える． また，温室ガス固定化や再利用化プロセスの開発においても，水溶液中のCO₂やH₂CO₃の状態は重要である． CO₂が水分子と反応してH₂CO₃となるのは知られているが，H₂CO₃の単離生成は実験的に難しく，低温ガス状態で初めてH₂CO₃作成に成功した [1]． その後液体窒素中に貯留したH₂CO₃の解離反応の反応定数決定により，H₂CO₃単体は安定物質である事が実験的に示された [2]． 量子化学計算によりH₂CO₃からCO₂への反応に水分子(n=1，2，3)が介在すると，反応活性化エネルギーが低下する事が明らかになり，水分子の存在がH₂CO₃の不安定性に影響を与えている事を実験と理論の両面により示された [3]．

CO₂，H₂CO₃および遷移状態に対する水分子を介した状態の理論計算は精力的に進められ，第一原理分子動力学(AIMD)計算によるH₂CO₃の解離状態の構造変化に関する解析が行われている [4, 5]． しかしながら，水分子ネットワーク中の様々なコンフォーメーションや反応機構の第一原理計算による探索には計算資源の観点で限界がある． 近年，様々な系で適用が報告されている機械学習ポテンシャル(MLP)は，エネルギーや原子間力をラベル値とした構造情報データセットを教師データとし教師アリ学習により生成する原子間ポテンシャルである [6]． 今回，DeepMDフレームワーク [7,8,9]を用いてMLPを作成し，水分子，炭酸，CO₂によって構成される水分子ネットワークに対する分子動力学(MD)計算によるコンフォーマの解析，Nudged Elastic Band (NEB) 法による遷移状態探索を行い，AIMDと比較して有効性を検討した．

2 Method

教師データ作成は第一原理計算プログラムパッケージ(VASP5.4.1) [10,11,12]を用い，交換相関ポテンシャルとしてPerdew-Burke-Ernzerhof(PBE)を，分散力補正としてDFT-D3 [13]を取り入れた． 水分子のみ，溶質分子のみ，水分子+溶質分子(遷移状態類似構造含む)，水分子(イオン)，水分子+溶質分子(水和構造)の5グループの教師データを，全体で4.05M frames作成した． ポテンシャル生成はDeepMD(2.1.0)とdpgenの組合せによるadaptive learningで生成した． MLPを用いたMD計算にはlammps(29 Sep 2021) [14]を使用した． 構造モデルの作成と結果の可視化はWinmostart(V11.2.2)とVESTA(3.5.7)を用いた．

3 Results and Discussion

Figure 1はH₂CO₃のtrans-transコンフォーメーションを含む構造に対する温度300KでのNVTアンサンブルでのトラジェクトリを示す． 青がMLP-MDの結果， 橙がAIMDの結果である． Figure1(a)はO-C-Oの結合角を示しており， 2つのO=C-O結合角は125度を中心に振動していることが分かる． またこの挙動はMLP-MDとAIMDでほぼ同様であった． Figure 1(b)は3つのO-C結合長を示し， O=C結合部は約1.25Åを中心とし， それ以外のO–C結合長はやや長い1.35Åを中心として振動している． このMDの初期構造でのMLP-MDとAIMDのエネルギー差は0.00475 eV/atom (0.109 kcal/mol/atoms)となった．

Figure 2にNEB法で求めた解離反応パスの動径分布関数(RDF)を示す． Figure 2(a)(b)に，それぞれ，C–O対とH-O対を示す． (0)から(9)の数字はFigure 2(c)に示す反応座標の番号である．Figure 2(b)の縦破線は通常の水分子の結合長(0.98Å)を示す． (2)から(6)の遷移状態では，C-O結合の遷移的な変化だけでなく，水分子内のZundel/Eigen型に類似した通常水分子よりも長いH–O結合の水和構造を確認した．

Figure 3にNEB法で計算したエネルギーバリア値を過去文献と含めてプロットした． 今回の計算は30H₂O分子で構成した系での計算である． AIMD計算のバリア値は9.76 kcal/molであり，MLPを用いたMD計算の値は14.59 kcal/molであった． このエネルギー差は4.83 kcal/mol(208 meV)である． 水分子ネットワークのH₂O分子数に対するバリア値の傾向を見ると，H₂O分子数の増大に対してバリア値が減少していく傾向は過去の量子化学計算 [3]およびAIMD [5]と同様の傾向であった．

Figure 1

Trajectories of the H₂CO₃ structure. (a): the O-C-O bond angles、 (b): the O-C bond lengths.

Figure 2

RDF of the reaction path of the dissociation process; (a):C-O pairs、 (b):H-O pairs、 (c):reaction coordinates.

Figure 3

The energy barrier as a function of number of H₂O.

References

• [1]   M. H. Moore, R. K. Khanna, Spectrochim. Acta A, 47, 255 (1991). doi:10.1016/0584-8539(91)80097-3
• [2]   W. Hage, K. R. Liedl, A. Hallbrucker, E. Mayer, Science, 279, 1332 (1998). doi:10.1126/science.279.5355.1332 PMID:9478889
• [3]   C. S. Tautermann, A. F. Voegele, T. Loerting, I. Kohl, A. Hallbrucker, E. Mayer, et al., Chemistry, 8, 66 (2002). PMID:11822465
• [4]   P. P. Kumar, A. G. Kalinichev, R. J. Kirkpatrick, J. Phys. Chem. B, 113, 794 (2009). doi:10.1021/jp809069g PMID:19108639
• [5]   M. Galib, G. Hanna, J. Phys. Chem. B, 115, 15024 (2011). doi:10.1021/jp207752m PMID:22053746
• [6]   J. Behler, M. Parrinello, Phys. Rev. Lett., 98, 146401 (2007). doi:10.1103/PhysRevLett.98.146401 PMID:17501293
• [7]   H. Wang, L. Zhang, J. Han, W. e, Comput. Phys. Commun., 228, 178 (2018). doi:10.1016/j.cpc.2018.03.016
• [8]   L. Zhang, J. Han, H. Wang, R. Car, W. e, Phys. Rev. Lett., 120, 143001 (2018). doi:10.1103/PhysRevLett.120.143001 PMID:29694129
• [9]   Y. Zhang, H. Wang, W. Chen, J. Zeng, L. Zhang, H. Wang, et al., Comput. Phys. Commun., 253, 107206 (2020). doi:10.1016/j.cpc.2020.107206
• [10]G. Kresse, J. Hafner, Phys. Rev. B, 47, 558 (1993). Phys. Rev. B, 49 14251 (1994).
• [11]   G. Kresse, J. Furthmüller, Comput. Mater. Sci., 6, 15 (1996). doi:10.1016/0927-0256(96)00008-0
• [12]   G. Kresse, J. Furthmüller, Phys. Rev. B Condens. Matter, 54, 11169 (1996). doi:10.1103/PhysRevB.54.11169 PMID:9984901
• [13]   S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys., 132, 154104 (2010). doi:10.1063/1.3382344 PMID:20423165
• [14]   A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, S. J. Plimpton, Comput. Phys. Commun., 271, 108171 (2022).