[计算化学]分子动力学笔记

本文为某计算机本科生的分子动力学学习笔记，在gpt4的辅助下，非体系化地整理相关生物、化学、统计力学知识。所有生成内容经过检查和调整，均直接代表本人观点。有科学性错误的话欢迎指教。

什么是分子动力学

定义： 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。

自己的理解：使用计算机模拟分子系统，探索其行为、微观和宏观性质。

两类方法

一类基于第一性原理，使用薛定谔方程计算所有原子的波函数随时间变化，能够提供非常精确的模拟结果，因为它不依赖于任何经验参数。然而，由于需要解决电子的量子力学问题，这种方法计算成本非常高，通常限制于较小的系统或者较短的时间尺度。

一类基于统计力学原理，使用经验势（如Lennard-Jones势或者更复杂的多体势）来描述原子或分子间的相互作用。在这种方法中，系统的宏观性质是通过对大量微观状态的统计来预测的。

本文主要学习第二种

基本框架

1. 系统初始化

• 选择模型系统：确定要模拟的分子或原子组成的系统。
• 定义初始位置：为系统中的所有粒子（原子或分子）指定初始位置，通常这是通过实验数据、随机分布或晶体结构来确定。
• 分配初始速度：根据温度给粒子分配初始速度。这些速度通常是从玻尔兹曼分布中随机选取的。

2. 势能函数的选择

• 选择适当的势能函数：这是描述粒子间相互作用的数学表达式。势能函数可以是简单的两体相互作用，如Lennard-Jones势，也可以是复杂的多体相互作用，如反应场势或者键合势。

3. 力的计算

• 计算力：在每一步模拟中，根据势能函数计算作用在每个粒子上的力。这些力是势能对粒子坐标的负梯度。

4. 积分方程的求解

• 时间积分：使用数值积分算法，如Verlet算法或其变种（比如速度Verlet算法），来更新粒子的位置和速度。这些算法允许从当前粒子的位置和速度计算出一小时间步后的位置和速度。

玻尔兹曼分布

玻尔兹曼分布描述了在热平衡状态下，系统中的粒子在不同能量水平上的分布情况。这一分布是由奥地利物理学家路德维希·玻尔兹曼在19世纪末提出的，它揭示了温度、粒子能量与粒子数之间的关系。

数学表达

玻尔兹曼分布可以用以下公式表示：

\[P(E) \propto e^{-\frac{E}{k_BT}} \]

其中：

• \(P(E)\) 是粒子具有能量\(E\)的概率密度。
• \(E\) 是粒子的能量。
• \(k_B\) 是玻尔兹曼常数，其值约为\(1.38 \times 10^{-23} \, \text{J/K}\)。
• \(T\) 是系统的绝对温度（单位是开尔文）。

在分子动力学中

在分子动力学模拟中，需要初始化系统，初始化的系统里的粒子的分布就需要满足玻尔兹曼分布。

自由能

在热力学和统计物理中，自由能是一个描述系统能量状态以及系统在给定条件下进行自发过程的能力的重要物理量。自由能有两种常见的形式：亥姆霍兹自由能（Helmholtz free energy）和吉布斯自由能（Gibbs free energy），它们分别适用于不同的热力学条件。

亥姆霍兹自由能 (A)

亥姆霍兹自由能（\(A\)）是在恒温、恒体积条件下系统的自由能，它定义为：

\[A = U - TS \]

其中：

• \(A\) 是亥姆霍兹自由能，
• \(U\) 是系统的内能，
• \(T\) 是绝对温度，
• \(S\) 是系统的熵。

亥姆霍兹自由能的减少表示在恒温恒体积下系统自发进行的过程。

吉布斯自由能 (G)

吉布斯自由能（\(G\)）是在恒温、恒压条件下系统的自由能，它定义为：

\[G = H - TS \]

或者根据\(H = U + PV\)（\(H\)是焓，\(P\)是压力，\(V\)是体积）也可以写作：

\[G = U + PV - TS \]

其中：

• \(G\) 是吉布斯自由能，
• \(H\) 是系统的焓，
• \(P\) 是压力，
• \(V\) 是体积，
• \(T\) 和 \(S\) 分别是系统的绝对温度和熵。

吉布斯自由能的减少表示在恒温恒压下系统自发进行的过程。

在分子动力学中

由克劳修斯定理知孤立系统的熵只增不减，也就是自由能只降不增。分子动力学认为，自由能降到最低时系统达到平衡态。

熵

这部分建议阅读张老师的博客：https://www.cnblogs.com/manuscript-of-nomad/p/17154461.html#_label2_1

系综理论

系综理论是统计物理学的一个核心概念，它提供了一种方法来处理和理解大量粒子系统的宏观性质。在经典物理学中，理想情况下，如果我们知道系统中所有粒子的初始位置和速度，我们可以使用牛顿的运动定律精确地预测系统的未来状态。然而，对于包含大量粒子的系统（如气体、液体或固体），这种方法是不切实际的。系综理论通过使用概率和统计方法来解决这个问题，它不是关注单个粒子的行为，而是研究系统的宏观性质。

系综的基本概念

系综是指一大群虚拟的、宏观上处于相同条件下的微观系统的集合。每个微观系统都是独立的，并且具有相同的宏观条件（如体积、能量、粒子数等）。通过研究这个集合的统计性质，我们可以得到系统的热力学性质。

主要类型的系综

系综理论中有三种基本的系综，它们分别对应于不同的物理条件：

1. 微正则系综（Microcanonical Ensemble）：适用于能量、体积和粒子数都恒定的孤立系统。在这个系综中，所有微观状态具有相同的能量，因此每个状态出现的概率相等。微正则系综的研究重点是熵和温度的关系。

2. 正则系综（Canonical Ensemble）：适用于体积、粒子数和温度恒定的封闭系统。这种系统与一个大热库处于热平衡。在正则系综中，不同微观状态具有不同的能量，状态出现的概率遵循玻尔兹曼分布。正则系综的研究重点是系统的自由能。

3. 巨正则系综（Grand Canonical Ensemble）：适用于温度、体积和化学势恒定的开放系统。这种系统可以与外界交换粒子。在巨正则系综中，除了能量分布，粒子数也是变量，系统的性质可以通过巨配分函数来描述。

在系综理论中，相空间是一个非常重要的概念，它为描述和分析系统的微观状态提供了一个强大的框架。相空间是一个抽象的多维空间，其中的每一点代表系统的一个可能的微观状态。对于一个由粒子组成的物理系统，每个粒子的位置和动量（或速度）共同定义了系统的微观状态。

相空间的定义

对于一个单一粒子的简单系统，相空间由粒子的位置和动量坐标构成。例如，如果我们考虑一个自由移动的单粒子在三维空间中，其相空间将由六个维度构成：三个位置坐标\((x, y, z)\)和三个与之对应的动量坐标\((p_x, p_y, p_z)\)。因此，这个单粒子的相空间是一个六维空间，每个点\((x, y, z, p_x, p_y, p_z)\)代表粒子的一个可能的微观状态。

对于一个包含\(N\)个粒子的系统，相空间的维度将是\(6N\)，因为每个粒子都有三个位置坐标和三个动量坐标。因此，相空间的每一点都代表了整个系统的一个可能状态，即所有粒子的位置和动量的一个特定组合。

能量平面

在经典力学中，一个系统的哈密顿量（Hamiltonian）\(H\)是描述系统能量的函数，它通常依赖于系统的广义坐标\(q_i\)和与之共轭的广义动量\(p_i\)。对于一个不受外力作用、孤立的系统，哈密顿量是守恒的，即不随时间变化。

在相空间中，系统的状态可以由一个点表示，该点的位置由广义坐标和广义动量的值确定。随着时间的推移，这个点在相空间中的轨迹（称为相轨迹）将遵循哈密顿方程：

\[\frac{dq_i}{dt} = \frac{\partial H}{\partial p_i}, \quad \frac{dp_i}{dt} = -\frac{\partial H}{\partial q_i} \]

这里，\(dq_i/dt\)和\(dp_i/dt\)分别是广义坐标\(q_i\)和广义动量\(p_i\)随时间的变化率。

对于一个哈密顿量保持不变的孤立系统，其相轨迹将被限制在相空间中能量为常数的曲面上，即满足\(H(q_i, p_i) = E\)的曲面，其中\(E\)是系统的总能量。这个曲面被称为能量曲面或能量壳层。

配分函数

配分函数在统计物理学的系综理论中扮演着核心角色，它是连接微观物理状态与宏观热力学性质的桥梁。配分函数的具体形式取决于所研究的系综类型（微正则、正则或巨正则系综），但在所有情况下，它都提供了一种计算系统宏观性质的方法。

正则系综的配分函数

在正则系综中，系统与一个大热库处于热平衡，能够交换能量，但粒子数和体积保持不变。这种情况下的配分函数被称为正则配分函数，通常用\(Z\)表示。

对于一个由\(N\)个粒子组成的系统，处于温度\(T\)下，其正则配分函数定义为：

\[Z = \sum_{i} e^{-\beta E_i} \]

其中，求和遍及系统所有可能的微观状态\(i\)，\(E_i\)是系统在第\(i\)个微观状态下的能量，\(\beta = \frac{1}{k_B T}\)，\(k_B\)是玻尔兹曼常数，\(T\)是绝对温度。

正则配分函数\(Z\)的重要性在于，它可以用来计算系统的各种热力学量。例如，系统的自由能\(F\)、内能\(U\)、熵\(S\)、热容\(C\)等都可以通过\(Z\)及其对温度的导数来表达：

• 亥姆霍兹自由能：\(F = -k_B T \ln Z\)
• 内能：\(U = -\frac{\partial \ln Z}{\partial \beta}\)
• 熵：\(S = -\left(\frac{\partial F}{\partial T}\right)_V\)
• 热容：\(C_V = \left(\frac{\partial U}{\partial T}\right)_V\)

巨正则系综的配分函数

在巨正则系综中，系统不仅能够与热库交换能量，还能与粒子库交换粒子。这种情况下的配分函数被称为巨正则配分函数，通常用\(\Xi\)表示。

巨正则配分函数定义为：

\[\Xi = \sum_{N=0}^{\infty} \sum_{i} e^{-\beta (E_{i,N} - \mu N)} \]

其中，第一个求和遍及所有可能的粒子数\(N\)，第二个求和遍及在给定\(N\)下所有可能的微观状态\(i\)，\(E_{i,N}\)是系统在有\(N\)个粒子且处于第\(i\)个微观状态下的能量，\(\mu\)是化学势。

巨正则配分函数\(\Xi\)同样可以用来计算系统的各种热力学量，如系统的压强、化学势、平均粒子数等。

在分子动力学模拟中

通常使用正则系综（恒温恒压）。给定相空间中的一个点，我们就能计算它的能量，进而找到能量平面的最低点，也就是稳定态。能量是动能+势能，因为我们都是恒温模拟，所以相空间内所有状态总动能是不变的，我们主要关心势能。计算给定状态的总势能的方法被称为力场。

力场

分子动力学模拟的力场定义了模拟中分子间和分子内部各原子之间的相互作用力。力场的选择和参数化对模拟结果的准确性和可靠性有着决定性影响。

力场的组成

力场通常由以下几个部分组成：

1. 键合相互作用：这些相互作用描述了原子间的化学键。它们通常包括键长（bond lengths）、键角（bond angles）和二面角（dihedral angles）的势能函数。这些势能函数通常采用简单的数学形式，如谐振子势（对于键长和键角）和周期势（对于二面角）。

2. 非键合相互作用：这部分描述了不通过化学键直接相连的原子之间的相互作用，主要包括范德华力（Van der Waals forces）和库仑力（Coulomb forces）。范德华力通常使用Lennard-Jones势来描述，而库仑力则描述了带电粒子之间的静电相互作用。

势能函数

力场通过势能函数（Potential Energy Function）来描述这些相互作用。总势能\(U\)是各种相互作用势能的总和：

\[U = U_{\text{bond}} + U_{\text{angle}} + U_{\text{dihedral}} + U_{\text{nonbonded}} \]

其中，\(U_{\text{bond}}\)、\(U_{\text{angle}}\)、\(U_{\text{dihedral}}\)分别代表键合相互作用中的键长、键角和二面角的势能，而\(U_{\text{nonbonded}}\)则包括所有非键合相互作用的势能。

力场的参数化

力场的参数化是指为势能函数中的各个项确定参数（如力常数、平衡键长、平衡键角等）的过程。这些参数通常通过实验数据和量子化学计算来确定。准确的参数化对于模拟的成功至关重要，因为它直接影响到模拟所预测的分子结构、动力学和热力学性质的准确性。

粗粒化力场

全原子力场计算系统内所有原子的相互作用，进而求出精确的势能。而粗粒化力场使用一些策略将几个临近的原子打包成一个大粒子，然后使用人工制造的力场来求势能。这个势能也是人工的，它并不代表自然界真实的势能，它为研究目的服务。

Upside-md方法

upside-md是jumper的博士毕业论文。jumper也就是Alphafold2的一作，AF2是目前最具革命性的AI模型之一，被认为可以准确预测蛋白质的三级结构，在结构生物、制药等领域有革命性的影响。

upside-md是用于研究蛋白质动力学的模型，它的核心思想是抓住蛋白的主要矛盾（主链结构），将侧链打包成有方向的珠子。整个蛋白被数学表示为主链坐标和侧链chi角的序列，大大减少了粒子数。在给定主链坐标的情况下，侧链chi角的所有状态构成一个相空间。

使用参数化的侧链-侧链相互作用、侧链-主链相互作用力场计算相空间内状态的势能、进而求改主链坐标配置的自由能。这个参数使用机器学习方法最大似然方法求得，使用电镜测定的真实蛋白质结构作为训练集，训练目标是使相空间内，真实状态存在的概率尽量大（也就是真实状态能量尽量小，虚假状态能量尽量高）。

认为自由能是降低的，因此自由能对主链原子坐标的偏导数，就是该原子受到的力。然后运行朗之万动力学（保证系统始终满足玻尔兹曼分布），迭代多次后收敛到稳定状态。

该模型的计算量非常小，可以使用单核CPU在几天内收敛一个蛋白。可以用于蛋白的从头折叠（只给氨基酸序列，预测蛋白结构），蛋白的动力学分析（如多个平衡态之间变化的分析），蛋白质结合分析（蛋白-蛋白互作用、蛋白-小分子结合。