分子集合体の弾性体モデルにもとづく粗視化剛性行列の力学的解釈

Abstract

We have been developing a method for coarse-graining the low-frequency vibration modes of molecular assemblies, which affords a numerical representation of the down-sized stiffness matrix. In this study, we present an analytical representation of the stiffness matrix based on the elastic-body modeling of molecular assemblies. Comparison between the numerical and analytical data allows the 13 parameters regarding the dimension and mechanical properties of the putative elastic body. The results for 57 molecular dimers with various hydrogen-bond multiplicity demonstrate that the obtained parameters were physically reasonable and well-reproduced the wavenumbers of normal-mode vibrations.

Translated Abstract

We have been developing a method for coarse-graining the low-frequency vibration modes of molecular assemblies, which affords a numerical representation of the down-sized stiffness matrix. In this study, we present an analytical representation of the stiffness matrix based on the elastic-body modeling of molecular assemblies. Comparison between the numerical and analytical data allows the 13 parameters regarding the dimension and mechanical properties of the putative elastic body. The results for 57 molecular dimers with various hydrogen-bond multiplicity demonstrate that the obtained parameters were physically reasonable and well-reproduced the wavenumbers of normal-mode vibrations.

1 はじめに

結晶やタンパク質，核酸などの分子集合体における構造–機能相関を計算化学的に扱うには，精度の高い分子間力パラメータが必須である．静的構造の評価においては分子間力のポテンシャル関数の深さが重要である一方，動的構造の評価にはその曲率，すなわち剛性定数が必要となる．室温付近において動的構造に大きく寄与するのはTHzオーダーの低振動数(LF)モードの振動であるが，このモードを含む基準振動計算を，全原子を量子化学レベルで考慮して行うには大きな計算機負荷が伴う．また，この振動数領域では並進や秤動など分子の相対配置の変位振動に分子のたわみ・ねじれなどの分子内変位が結合するため，単純な剛体近似によるモデル化は難しい．我々はこれまでに，分子集合体の基準振動を粗視化してHesse行列の次元を削減するスキームを構築してきた [1,2,3,4]．この方法では分子二量体の剛性定数が最小で12×12の行列要素で表現されるが，要素間には一定の関係があるため全てが独立ではない．本研究では個々の要素の物理的な意味を明確にし，さらに少数のパラメータで剛性行列を近似することを目的として，二量体の力学的性質を古典的な弾性体力学の用語で解釈することを試みた．

2 方法

分子集合体を連成調和振動子とみなすと，その運動方程式は質量行列Mと剛性行列Kからなる質量加重Hesse行列の固有値問題に帰せられる．分子集合体の構成原子数をNとすると行列の次元は3Nであり，非直線分子では6個の0固有値と3N−6個の非0固有値を与える．ここで原子の変位を各構成分子の3個の純並進，3個の純回転，α個の分子内振動モードに変換する粗視化行列Bを考えると，式1,2のΓ⁻¹とΦはそれぞれ粗視化された空間における質量と剛性の表現行列となる [1, 2]．   

${\Gamma }^{-1}=B^{T}MB$(1)

  

$\Phi =B^{T}KB$(2)

運動方程式を表す固有値方程式は，固有ベクトルU，固有値ω²を要素とする対角行列Ω²もちいて以下のように書ける．   

${\Gamma }^{1/2}\Phi {\Gamma }^{1/2}U=U{\Omega }^2$(3)

基準振動解析の結果(U, Ω²)があれば，あらわにKを求めなくとも式3の逆変換によりΦ が計算できる．   

$\Phi ={\Gamma }^{-1/2}U{\Omega }^2{U}^T{\Gamma }^{-1/2}$(4)

次にΦの力学的な表現式について考える．分子二量体のモデルとして，二個の剛体(unit I, II)が結合した直方体の弾性体を仮定した(Figure 1)．

Figure 1.

 Elastic-body model for the molecular dimer.

弾性体の中心を原点におき，長さL, W, Hの辺がx, y, z軸方向に一致するように座標系(global系)をとった．この座標系でunit I, IIの重心をそれぞれ(a₁, b₁, c₁), (a₂, b₂, c₂)とし，これらを原点とする座標系(local-I系, -II系)を定義した．弾性体の頂点1–4をlocal-I系で，頂点5–8をlocal-II系で表し，式1, 2と同様の考え方で剛体の純並進・純回転にともなう変位を24×12行列B_Eで表した．これらの変位を弾性体のひずみに変換する9×24行列をDとすることにより，粗視化剛性行列と等価な行列の表式として，   

${\Phi }_{\text{E}}=V{B}_{\text{E}}{}^T{D}^{T}ED{B}_{\text{E}}$(5)

を導出した．ここでVは弾性体の体積(= LWH)，Eは弾性定数などからつくられる9×9行列である(式6, 7)．   

$E=\left(\begin{array}{ccc}\lambda & & O\\ & \mu & \\ O& & \rho \end{array}\right)$(6)

  

$\{\begin{array}{c}\lambda =\left(\begin{array}{ccc}{\lambda }_x+2{\mu }_{xx}& {\lambda }_x& {\lambda }_x\\ {\lambda }_y& {\lambda }_y+2{\mu }_{yy}& {\lambda }_y\\ {\lambda }_z& {\lambda }_z& {\lambda }_z+2{\mu }_{zz}\end{array}\right),\\ \mu =\left(\begin{array}{ccc}{\mu }_{yz}& 0& 0\\ 0& {\mu }_{zx}& 0\\ 0& 0& {\mu }_{xy}\end{array}\right),\rho =\left(\begin{array}{ccc}{\rho }_x& 0& 0\\ 0& {\rho }_y& 0\\ 0& 0& {\rho }_z\end{array}\right)\end{array}$(7a–c)

ここでλ_s, μ_st (s, t = x, y, z)は異方性体に拡張したラメの定数で，ヤング率(E_x)，横弾性係数(G_xy, G_zx)，ポアソン比(ν)と以下の関係にある．   

$\{\begin{array}{c}{\lambda }_s=\frac{\nu E_s}{(1+\nu )(1-2\nu )},\\ {\mu }_{ss}=\frac{E_s}{2(1+\nu )},{\mu }_{st}=G_{st}\end{array}$(8a–c)

またρ_s (s = x, y, z)は断面二次極モーメント(I_p)，断面二次モーメント(I_y, I_z)などの情報を含む．   

${\rho }_x=\frac{G_{yz}{I}_p}{LV},{\rho }_y=\frac{E_x{I}_y}{LV},{\rho }_z=\frac{E_x{I}_z}{LV}$(9a–c)

  

$\{\begin{array}{c}{I}_p=I_y+I_z\\ I_y={\int }_0^W{\int }_{-H/2}^{H/2}{z}^2\text{d}z\text{d}y=\frac{H^{3}W}{12}\\ I_z={\int }_0^H{\int }_{-W/2}^{W/2}{y}^2\text{d}y\text{d}z=\frac{H{W}^3}{12}\end{array}$(10a–c)

D行列を全て書き尽くすには紙面が足りないが，例えばx軸方向の引張ひずみについての行はiを頂点の番号として式11で与えられる．   

$\begin{array}{c}{\left(}_D=\left(\begin{array}{cc}{d}_1\cdots d_4& d_5\cdots d_8\end{array}\right)\\ d_i=\{\begin{array}{c}{\left(}^4\left(\begin{array}{ccc}1& 0& 0\end{array}\right),i=\{1,2,3,4\}\\ {\left(}^4\left(\begin{array}{ccc}-1& 0& 0\end{array}\right),i=\{5,6,7,8\}\end{array}\end{array}$(11)

Φ_Eの要素，例えば並進に関係するブロック{Φ_E }_TTは式12のようになる．   

${\left\{{\Phi }_{\text{E}}\right\}}_{\text{TT}}=\frac{V}{L^2}\left(\begin{array}{ccc}{\lambda }_x+2{\mu }_{xx}& 0& 0\\ 0& {\mu }_{xy}& 0\\ 0& 0& {\mu }_{zx}\end{array}\right)$(12)

水素結合により会合した57種類の分子二量体をHF/6-311G**条件で構造最適化し，次いで基準振動解析を行った．計算にはGaussian09および16を用いた [5, 6]．得られた振動数と原子変位のデータから，式(4)に基づいてΦ行列を計算した†．この粗視化処理で得られたΦを，式(5)のΦ_Eと比較することにより，a_1,b₁, c₁, a₂, b₂, c₂, E_x, E_y, E_z, ν, L, W, Hの13個のパラメータを算出した．ただしLだけは係数比較からは直接求められなかったため，水素結合に関与するヘテロ原子間の距離の平均値を採用した．

3 結果と考察

N-メチルアデニン–N-メチルチミン対の結果を例として示す．この系のLFモード基準振動は，純並進・純回転の6次元に加えて分子内振動モードを2個考慮(α = 2)した14次元の粗視化空間でよく再現することができた．弾性体モデルへの回帰によって得られたパラメータは，(a₁, b₁, c₁) = (3.684, −0.501, −0.004) Å，(a₂, b₂, c₂) = (−2.882, −0.704, 0.000) Å，E_x = 115.2 GPa, G_xy = 52.0 GPa, G_zx = 11.1 GPa, ν = 0.305, (L, W, H) = (3.05, 4.81, 1.21) Åとなった．local-I, -II系の原点間の距離は6.569 Åであり，構成分子の重心間の距離6.564 Åとよく一致した．global系の原点は水素結合部位の中央付近となった．Figure 2に，弾性体を直方体で，unit I, IIを各構成分子の慣性テンソルと等価な楕円体で示す．

Figure 2.

 Visualization of the result of regression to the elastic-body model of the molecular dimer.

これらの弾性体パラメータから剛性行列を再構成し，式(3)に基づいて固有振動数を計算したところ，Gaussianによる基準振動計算の結果とよい一致を示した(Figure 3)．これは，水素結合による会合体の力学的性質が古典的な弾性体に基づいて解釈できることを表している．

Figure 3.

 Correlation of wavenumbers of LF-mode vibrations between the original normal-mode calculation and results of the coarse-graining analyses.

Figure 2の座標系設定下では，力の対称性のためΦ行列の並進・回転に関係する12×12 = 144個の行列要素のうち6×6 = 36個が独立である．本研究の弾性体近似によって独立なパラメータは13個にまで圧縮され，物理的意味も明確になった．

他の分子二量体についても同様の解析を行い，その多くは基準振動の固有振動数が弾性体モデルによってよく再現された．また，水素結合の多重度が増すにしたがって弾性体のy軸方向の長さにあたるWは大きくなり，それに伴って(V/L²)E_x（力の定数に相当）も増加した(Figure 4)．一方，2-アミノピリジン二量体など数種類の系についてはポアソン比が許容範囲外になるなどして，弾性体モデルへの回帰自体ができなかった．

Figure 4.

 Correlation between the W value and Yang modulus (E_x). The plots of molecular dimers with the same hydrogen bond multiplicity are grouped by labeled circles. The molecular structures with a cuboid are the representatives of each group.

今後は弾性体モデルの適用範囲を明確にするとともに，多量体や周期系 [3, 4]に応用するための拡張を試みる．分子集合体を剛体と弾性体の結合で近似することにより，トランススケールの分子レベル有限要素法シミュレーションが実行可能になると期待される．

†The program employed is not in public as it is under development. Those who are interested in co-developing are welcome to contact the corresponding author.

参考文献

• [1]   H. Houjou, J. Chem. Theory Comput., 5, 1814 (2009). doi:10.1021/ct900169f PMID:26610005
• [2]   H. Houjou, J. Chem. Phys., 135, 154111 (2011). doi:10.1063/1.3652102 PMID:22029301
• [3]   H. Houjou, Y. Wang, S. Okamura, J. Comput. Chem. Jpn., 20, 147 (2021). doi:10.2477/jccj.2022-0004
• [4]   H. Houjou, M. Seshimo, J. Math. Chem., 60, 613 (2022). doi:10.1007/s10910-022-01327-w
• [5]   M. J. Firsch, et al., Gaussian09 Rev. D.01 Gaussian, Inc., Wallingford CT, 2013.
• [6]   M. J. Firsch, et al., Gaussian16 Rev. C.01 Gaussian, Inc., Wallingford CT, 2019.