关注 $\mathrm{M}\mathrm{r}.\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l},$ $更$$多$$精$$彩$！

$主要专栏内容包括：$

$\text{ †《LAMMPS小技巧》：}$ 主要介绍采用分子动力学（$Lammps$）模拟相关安装教程、原理以及模拟小技巧（难度：$★$）
$\text{ ††《LAMMPS实例教程—In文件详解》：}$ 主要介绍采用分子动力学（$Lammps$）模拟相关物理过程模拟。（包含：热导率计算、定压比热容计算，难度：$★$$★$$★$）
$\text{ †††《Lammps编程技巧及后处理程序技巧》：}$ 主要介绍针对分子模拟的动力学过程（轨迹文件）进行后相关的处理分析（需要一定编程能力。难度：$★$$★$$★$$★$$★$）。
$\text{ ††††《分子动力学后处理集成函数—Matlab》：}$ 主要介绍针对后处理过程中指定函数，进行包装，方便使用者直接调用（需要一定编程能力，难度：$★$$★$$★$$★$）。
$\text{ †††††《SCI论文绘图—Python绘图常用模板及技巧》：}$ 主要介绍针对处理后的数据可视化，并提供对应的绘图模板（需要一定编程能力，难度：$★$$★$$★$$★$）。

漫谈自由能—未整理好请先不要看

• 一、理论篇——自由能要算什么？
• • 配分函数的 $Boltzmann$ 关系定义为 $Helmholtz$ 自由能
  • 配分函数的 $Boltzmann$ 关系定义为 $Gibbs$ 自由能
  • • 对于 $NVT$ 的配分函数
    • QNVT=1N!h3N∬eH(r⃗,p⃗)kBTdNr⃗dNp⃗Q_{NVT}=\frac{1}{N!h^{3N}}\iint e^\frac{{H(\vec{r},\vec{p})}}{k_BT}d^{N}\vec{r}d^{N}\vec{p}QNVT​=N!h3N1​∬ekB​TH(r,p​)​dNrdNp​
    • 对于 $NPT$ 的配分函数
    • QNPT=1N!h3N∭eH(r⃗,p⃗)pVkBTdNr⃗dNp⃗dVQ_{NPT}=\frac{1}{N!h^{3N}}\iiint e^\frac{{H(\vec{r},\vec{p})pV}}{k_BT}d^{N}\vec{r}d^{N}\vec{p}dVQNPT​=N!h3N1​∭ekB​TH(r,p​)pV​dNrdNp​dV
• 二、实践篇——自由能要如何计算？
• • 微扰动法
  • 热力学积分法
  • 平均力势法
• 三、总结

一、理论篇——自由能要算什么？

配分函数的 $Boltzmann$ 关系定义为 $Helmholtz$ 自由能

F=−kBTInQNVTF=-k_{B}TInQ_{NVT}F=−kB​TInQNVT​（$NVT$）

配分函数的 $Boltzmann$ 关系定义为 $Gibbs$ 自由能

F=−kBTInQNPTF=-k_{B}TInQ_{NPT}F=−kB​TInQNPT​（$NPT$）

如果需要计算能量，我们只需要计算当前构型全部粒子的结构就可以获得能量。
那么对于自由能，我们需要知道熵，S=kBInΩS=k_{B}In \OmegaS=kB​InΩ，$\Omega$ ：就是有多少种微观状态分布。那么怎么获得这个微观状态数呢？以$NVT$为例子，即 F=−kBTInQNVTF=-k_{B}TInQ_{NVT}F=−kB​TInQNVT​。

如下体系在 $NVT$ 系综下进行动力学模拟。$t$ 时间后, 输出每一帧的轨迹，也就是每一帧都对应的一个构型。

那么对于宏观上具有相同 $NVT$ 的结构，具有多少微观构型呢？即宏观上相同，但是微观上不同，这样的体系有多少个呢？
那么在这么一个连续分布、粒子数比较多的情况，相同的$NVT$ 有多少个呢？可以理解每输出一个就是一个不同微观状态的 $NVT$ 那么无限长时间能都遍历全部的 $NVT$ 微观状态呢？
根据庞加莱回归（see below）理论上是可以以一个很小的误差回归到一个状态，但是这个时间可能是个天文数字。我们无法穷举。

对于 $NVT$ 的配分函数

QNVT=1N!h3N∬eH(r⃗,p⃗)kBTdNr⃗dNp⃗Q_{NVT}=\frac{1}{N!h^{3N}}\iint e^\frac{{H(\vec{r},\vec{p})}}{k_BT}d^{N}\vec{r}d^{N}\vec{p}QNVT​=N!h3N1​∬ekB​TH(r,p​)​dNrdNp​

显然，对于 dNr⃗dNp⃗d^{N}\vec{r}d^{N}\vec{p}dNrdNp​ 所有可能的 r⃗\vec{r}r 和 p⃗\vec{p}p​ 相空间的积分，我们是无法直接得到的。

对于 $NPT$ 的配分函数

QNPT=1N!h3N∭eH(r⃗,p⃗)pVkBTdNr⃗dNp⃗dVQ_{NPT}=\frac{1}{N!h^{3N}}\iiint e^\frac{{H(\vec{r},\vec{p})pV}}{k_BT}d^{N}\vec{r}d^{N}\vec{p}dVQNPT​=N!h3N1​∭ekB​TH(r,p​)pV​dNrdNp​dV

那么如何化简推导呢？

H(r⃗,p⃗)=12m∑p⃗2+U(r⃗)H(\vec{r},\vec{p})=\frac{1}{2m}\sum \vec{p}^{2}+U(\vec{r})H(r,p​)=2m1​∑p​2+U(r)

那么 可以拆开写为动能的贡献部分：12m∑p⃗2\frac{1}{2m}\sum \vec{p}^{2}2m1​∑p​2；以及势能的贡献部分：U(r⃗)U(\vec{r})U(r)

QNVT=1N!h3N∬eH(r⃗,p⃗)kBTdNr⃗dNp⃗Q_{NVT}=\frac{1}{N!h^{3N}}\iint e^\frac{{H(\vec{r},\vec{p})}}{k_BT}d^{N}\vec{r}d^{N}\vec{p}QNVT​=N!h3N1​∬ekB​TH(r,p​)​dNrdNp​

上式就可拆开为:

QNVT=1N!h3N∬e12m∑p⃗2+U(r⃗)kBTdNr⃗dNp⃗Q_{NVT}=\frac{1}{N!h^{3N}}\iint e^\frac{{\frac{1}{2m}\sum \vec{p}^{2}+U(\vec{r})}}{k_BT}d^{N}\vec{r}d^{N}\vec{p}QNVT​=N!h3N1​∬ekB​T2m1​∑p​2+U(r)​dNrdNp​

动能部分可以直接积分计算：

Qk=1N!h3N∬e12m∑p⃗2kBTdNp⃗=Q_{k}=\frac{1}{N!h^{3N}}\iint e^\frac{{\frac{1}{2m}\sum \vec{p}^{2}}}{k_BT}d^{N}\vec{p}=Qk​=N!h3N1​∬ekB​T2m1​∑p​2​dNp​= 1N!Λ3N\frac{1}{N!\Lambda^{3N}}N!Λ3N1​, 其中，Λ=frach2πmkBT\Lambda=frac{h}{\sqrt{2\pi m k_BT} }Λ=frach2πmkB​T​

因此，对于$NVT$来说，可以写为：
QNVT=1N!Λ3N∬eU(r⃗)kBTdNr⃗Q_{NVT}=\frac{1}{N!\Lambda^{3N}}\iint e^\frac{{U(\vec{r})}}{k_BT}d^{N}\vec{r}QNVT​=N!Λ3N1​∬ekB​TU(r)​dNr

因此，对于$NPT$来说，可以写为：
QNPT=1N!Λ3N∬e−U(p⃗)kBTe−pvkBTdNr⃗dVQ_{NPT}=\frac{1}{N!\Lambda^{3N}}\iint e^\frac{{-U(\vec{p})}}{k_BT}e^\frac{{-pv}}{k_BT}d^{N}\vec{r}dVQNPT​=N!Λ3N1​∬ekB​T−U(p​)​ekB​T−pv​dNrdV

二、实践篇——自由能要如何计算？

微扰动法

热力学积分法

平均力势法

三、总结