空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

基于热力学积分(TI)的自由能计算

notebook

分子动力学

科学计算

DeePMD

Molecular Dynamics

bohrium

NBHub

Tutorial

LAMMPS

力场

notebook分子动力学科学计算DeePMDMolecular DynamicsbohriumNBHub TutorialLAMMPS力场

yfb222333@gmail.com

发布于 2023-12-22

推荐镜像 :Basic Image:ubuntu22.04-mojo0.3-notebook

推荐机型 :c2_m4_cpu

基于热力学积分的自由能计算

自由能计算

分子动力学模拟

热力学积分

基于热力学积分的自由能计算

这篇notebook将带介绍一下计算自由能的一个经典方法——热力学积分方法（Thermodynamics Integration方法）：并提供一个谐振子模型，计算谐振子变化到指定势函数时的自由能的变化。

bohrium上，已经有比较多notebook阐述了自由能计算方方面面，下面是节选的一些：

自由能计算

关于自由能，你想了解的一切（上）
**一句话导读：**对自由能计算FEP, BAR/MBAR, 热力学积分(TI),伞形采样 (Umbrella sampling)基本介绍
DP-GEN+FEP-TI计算自由能
**一句话导读：**介绍了一个工作流，训练势函数，计算MD中自由能的变化，并可视化
非平衡采样的自由能计算
**一句话导读：**对非平衡采样的自由能计算这一主题进行了非常系统的介绍，从原理到demo一应俱全
基于Torchmd实现伞型采样计算自由能
**一句话导读：**基于TorchMD包，介绍了伞形采样，并对丙氨酸二肽Phi/Psi旋转的PMF的一个窗口进行了计算。

分子动力学模拟

AI4S时代玩转分子动力学？看看你能到第几层！
**一句话导读：**介绍了一下分子动力学模拟基本概念，经典谐振子体系，保辛（相空间体积）的Velocity Verlet算法。

代码

文本

热力学积分

这一原理和形式足够清晰。但具体到计算时，需要额外考虑的需要有：数值积分计算，对误差（热力学误差，积分误差）估计，合理积分路径。

考虑两个有着不同哈密顿量的系统， $H_{0}$ 与 $H_{1}$ 我们想要计算这两个系统之间得自由能差，考虑一系列中间状态，其哈密顿量 $H (λ)$ 。随着 $λ = 0$ 逐渐变化到 $λ = 1$ ，系统由 $H_{0}$ 变化到 $H_{1}$ .

这两个系统之间的自由能差 $Δ A$ 满足， $Δ A = 0 \int 1 \frac{\partial A ( λ )}{\partial λ} d λ$

= $0 \int 1 - β^{- 1} \frac{\partial l n Q ( λ )}{\partial λ} d λ$

= $0 \int 1 - \frac{β ^{- 1}}{Q ( λ )} \frac{\partial Q ( λ )}{\partial λ} d λ$

= $0 \int 1 - \frac{β ^{- 1}}{Q ( λ )} \frac{\partial \int e ^{- β H (x, λ)} d x}{\partial λ} d λ$

= $0 \int 1 \frac{\int \frac{\partial H ( λ )}{\partial λ} e ^{- β H (x, λ)} d x}{Q ( λ )} d λ$

= $0 \int 1 ⟨ \frac{\partial H ( λ )}{\partial λ} ⟩_{λ} d λ$

其中 $Q (λ)$ 为不同 $λ$ 下，系统的配分函数。尖括号对 $<>$ 代表系综平均。

不同 $λ$ 对应系统不同的势函数。跑MD，实际上就是在这个势函数下，采样我们关心的物理量。也就是用MD轨迹上的点，代表系综平均值。

具体计算时，想要求得自由能差值，需要知道，中间这一系列状态的这一物理量的MD模拟的系综平均值，取路径上一系列积分节点，求数值积分即可。

一个简单的做法是让H(\lambda) 成为 $H_{0}$ 与 $H_{1}$ 的线性组合，即

也可以有一些更复杂的路线了。比如引入一些中间状态。

$H (λ) = (1 - λ) H_{0} + λ H_{1}$

这样这个 $\frac{\partial H ( λ )}{\partial λ}$ 就有相对简单的形式，在MD模拟过程中，可以简单、直接地算出来 $\frac{\partial H ( λ )}{\partial λ} = H_{1} - H_{0}$

MD软件实现时，运动轨迹的每一帧，得到每个原子位置后，再多算两个 $H_{0}$ 与 $H_{1}$ 两个势函数下系统的能量，并作差。

因此，自由能差的计算成为 \lambda 空间的积分。用最简单的积分法，梯形积分，它变成

$Δ A = i \sum \frac{1}{2} (⟨ H_{1} - H_{0} ⟩_{λ_{i}} + ⟨ H_{1} - H_{0} ⟩_{λ_{i + 1}}) Δ λ$

当然也可以由其他形式的更准确的数值积分方法。比如经典Runge-Kutta积分方法。

起始点一般可以选用爱因斯坦固体和理想气体，这两个状态自由能我们是可以解析求解。这样我们就可以得到目标系统自由能的具体数值。

代码

文本

双击即可修改

代码

文本

双击即可修改

代码

文本

[57]

!pip install numba

!pip install seaborn

!pip install scipy

Requirement already satisfied: numba in /opt/conda/lib/python3.8/site-packages (0.58.1)
Requirement already satisfied: llvmlite<0.42,>=0.41.0dev0 in /opt/conda/lib/python3.8/site-packages (from numba) (0.41.1)
Requirement already satisfied: numpy<1.27,>=1.22 in /opt/conda/lib/python3.8/site-packages (from numba) (1.24.4)
Requirement already satisfied: importlib-metadata in /opt/conda/lib/python3.8/site-packages (from numba) (6.8.0)
Requirement already satisfied: zipp>=0.5 in /opt/conda/lib/python3.8/site-packages (from importlib-metadata->numba) (3.17.0)
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv
Requirement already satisfied: seaborn in /opt/conda/lib/python3.8/site-packages (0.13.0)
Requirement already satisfied: numpy!=1.24.0,>=1.20 in /opt/conda/lib/python3.8/site-packages (from seaborn) (1.24.4)
Requirement already satisfied: pandas>=1.2 in /opt/conda/lib/python3.8/site-packages (from seaborn) (2.0.3)
Requirement already satisfied: matplotlib!=3.6.1,>=3.3 in /opt/conda/lib/python3.8/site-packages (from seaborn) (3.7.3)
Requirement already satisfied: contourpy>=1.0.1 in /opt/conda/lib/python3.8/site-packages (from matplotlib!=3.6.1,>=3.3->seaborn) (1.1.1)
Requirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.8/site-packages (from matplotlib!=3.6.1,>=3.3->seaborn) (0.11.0)
Requirement already satisfied: fonttools>=4.22.0 in /opt/conda/lib/python3.8/site-packages (from matplotlib!=3.6.1,>=3.3->seaborn) (4.42.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /opt/conda/lib/python3.8/site-packages (from matplotlib!=3.6.1,>=3.3->seaborn) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /opt/conda/lib/python3.8/site-packages (from matplotlib!=3.6.1,>=3.3->seaborn) (23.0)
Requirement already satisfied: pillow>=6.2.0 in /opt/conda/lib/python3.8/site-packages (from matplotlib!=3.6.1,>=3.3->seaborn) (10.0.1)
Requirement already satisfied: pyparsing>=2.3.1 in /opt/conda/lib/python3.8/site-packages (from matplotlib!=3.6.1,>=3.3->seaborn) (3.1.1)
Requirement already satisfied: python-dateutil>=2.7 in /opt/conda/lib/python3.8/site-packages (from matplotlib!=3.6.1,>=3.3->seaborn) (2.8.2)
Requirement already satisfied: importlib-resources>=3.2.0 in /opt/conda/lib/python3.8/site-packages (from matplotlib!=3.6.1,>=3.3->seaborn) (6.1.0)
Requirement already satisfied: pytz>=2020.1 in /opt/conda/lib/python3.8/site-packages (from pandas>=1.2->seaborn) (2023.3.post1)
Requirement already satisfied: tzdata>=2022.1 in /opt/conda/lib/python3.8/site-packages (from pandas>=1.2->seaborn) (2023.3)
Requirement already satisfied: zipp>=3.1.0 in /opt/conda/lib/python3.8/site-packages (from importlib-resources>=3.2.0->matplotlib!=3.6.1,>=3.3->seaborn) (3.17.0)
Requirement already satisfied: six>=1.5 in /opt/conda/lib/python3.8/site-packages (from python-dateutil>=2.7->matplotlib!=3.6.1,>=3.3->seaborn) (1.16.0)
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv
Collecting scipy
  Downloading scipy-1.10.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (34.5 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 34.5/34.5 MB 5.3 MB/s eta 0:00:0000:0100:01
Requirement already satisfied: numpy<1.27.0,>=1.19.5 in /opt/conda/lib/python3.8/site-packages (from scipy) (1.24.4)
Installing collected packages: scipy
Successfully installed scipy-1.10.1
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv

代码

文本