本文主要是介绍使用Metropolis蒙特卡洛方法的原子模拟，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
文章目录
- 1.蒙特卡罗方法的目标
- 2.热力学系综
- 3.连续体系
- 4.Metropolis算法
- 1.Metropolis算法介绍
- 2.Metropolis算法思路
- 5.原子体系的蒙特卡洛算法
- 1.算法的基本思想
- 2.算法的实现过程
1.蒙特卡罗方法的目标
- 蒙特卡罗方法可以做什么？
提供材料的热力学信息；
评估整体的平均值（能量、压力等）。 - 蒙特卡罗方法不能做什么？
提供材料的动力学信息（例如：扩散常数）
2.热力学系综
- 强度性质 广延性质
强度性质(intensive property)是不随物质多少或系统大小而改变的物理性质。
广延性质(extensive property)是一种物理性质，系统中此性质的量可由组成此系统所有子系统中对应性质的量相加而得。 - 热力学系综定义了哪些热力学量是受限制的，以及哪些热力学量是自由变化的。
- 微正则系综
在微正则系综（NVE）中，系统的所有可能状态都具有相同的能量 E，因此都具有相同的概率。
3.连续体系
- 对于原子或分子系统，微观状态的数量几乎是无限的，因此将总和替换为一组积分：
ρ ( r N ) = e − U ( r N ) / k B T Z N V T Z N V T = ∫ e − U ( r N ) / k B T d r N \rho(r^{N})=\frac{e^{-U(r^{N})}/k_{B}T}{Z_{NVT}}\quad\quad\quad Z_{NVT}=\int e^{-U(r^{N})/k_{B}T}dr^{N} ρ(rN)=ZNVTe−U(rN)/kBTZNVT=∫e−U(rN)/kBTdrN
例如：势能U的平均值:
< U > = 1 Z N V T ∫ e − U ( r N ) / k B T U ( r N ) d r N <U>=\frac{1}{Z_{NVT}}\int e^{-U(r^{N})/k_{B}T}U(r^{N})dr^{N} <U>=ZNVT1∫e−U(rN)/kBTU(rN)drN
为了评估这个物理量，应该先计算所有状态的能量并评估相应的概率。然而，大多数状态具有很高能量，并因此极不可能发生(这是一种浪费时间的做法)。
为了有效地评估该值，应该只对那些合理的少数状态进行求和，并排除其他状态。
4.Metropolis算法
1.Metropolis算法介绍
< U > = 1 Z N V T ∫ e − U ( r N ) / k B T U ( r N ) d r N <U>=\frac{1}{Z_{NVT}}\int e^{-U(r^{N})/k_{B}T}U(r^{N})dr^{N} <U>=ZNVT1∫e−U(rN)/kBTU(rN)drN
-
当我们需要求解一个 3N 维的积分时，
可以使用蒙特卡罗方法，通过在一个空间内进行随机抽样，该空间包含积分定义区域的体积，统计哪些点在该体积内，哪些点不在该体积内来实现。
ρ ( r N ) = e − U ( r N ) / k B T Z N V T \rho(r^{N})=\frac{e^{-U(r^{N})}/k_{B}T}{Z_{NVT}} ρ(rN)=ZNVTe−U(rN)/kBT -
Metropolis算法（1953年）是一种对构型空间进行采样的方法，其采样方式是以麦克斯韦分布给出的概率ρ "访问 "给定状态。
该方法的思路：在具有合理概率的状态集合中计算，得出该物理量的平均值。 -
Metropolis算法可以得到能量分布曲线，并提供具有合理概率的构型列表。得出的构型列表称为轨迹。
2.Metropolis算法思路
-
Metropolis算法的方法思路：
1.从给定的随机构型开始
2.进行一次试探性移动以获得新的构型
3.将新构型与之前构型进行比较，根据概率来判断是否将新构型添加到轨迹中。
ρ β ρ α = e − E β / k β T Q Q e − E α / k β T = e − ( E β − E α ) / k B T \frac{\rho_{\beta}}{\rho_{\alpha}}=\frac{e^{-E_{\beta}/k_{\beta}T}}{Q}\frac{Q}{e^{-E_{\alpha}/k_{\beta}T}}=e^{-(E_{\beta}-E_{\alpha})/k_{B}T} ραρβ=Qe−Eβ/kβTe−Eα/kβTQ=e−(Eβ−Eα)/kBT -
在轨迹上添加新构型的概率：
1.基于麦克斯韦分布的构型，其试探性移动（从配置i到配置i+1）将根据以下规则被接受：
当 Δ E i , i + 1 ≤ 0 \Delta E_{i,i+1}\le0 ΔEi,i+1≤0时，由于概率 e − Δ E i , i + 1 / k β T ≥ 1 e^{-\Delta E_{i,i+1}/k_{\beta}T}\ge1 e−ΔEi,i+1/kβT≥1而被接受；
当 Δ E i , i + 1 > 0 \Delta E_{i,i+1}>0 ΔEi,i+1>0时，以概率 e − Δ E i , i + 1 / k β T e^{-\Delta E_{i,i+1}/k_{\beta}T} e−ΔEi,i+1/kβT接受移动。
(可以通过生成介于0和1之间的随机数，并将其与设定的概率进行比较，来决定是否接受给定概率下的新构型。)
2.通过多次重复此过程，可以生成一系列构型 {i = 1…N}，其特征为能量 Ei，这些构型总体上具有合理的概率。 -
Metropolis算法的计算过程
1.初始化：从构型 i = 1 i=1 i=1开始，能量为 E i E_i Ei
2.进行一次随机试探性移动到构型 i + 1 i+1 i+1，能量为 E i + 1 E_{i+1} Ei+1，并计算能量差 Δ E \Delta E ΔE
如果 Δ E ≤ 0 \Delta E \le 0 ΔE≤0，则接受移动，将新构型添加到轨迹中。
如果 Δ E > 0 \Delta E >0 ΔE>0，则生成介于0和1之间的随机数 r r r。如果 r ≤ e x p ( − Δ E / k T ) r≤exp(-ΔE/kT) r≤exp(−ΔE/kT)，则接受移动；否则拒绝。
在高温下，可接受与当前状态能量差较大的新状态；在低温下，只接受与当前状态能量差较小的新状态。
3.如果移动被接受，新构型成为状态 i + 1 i+1 i+1。
如果移动被拒绝，构型 i + 1 i+1 i+1保持不变。每次试探都应该向轨迹中添加一个构型，尽管在这种情况下，它与前一个相同。
4. i = i + 1 i=i+1 i=i+1 ，转到步骤2。 -
蒙特卡洛方法的缺点
1.Metropolis蒙特卡罗方法是评估系统平均热力学的一种高效方式，然而，它并不提供系统动力学的任何信息。
2.没有明确的时间概念，“轨迹”不一定是系统自发遵循的轨迹。
然而，蒙特卡罗方法比分子动力学方法更通用，因为它可以应用于更广泛的系统。
5.原子体系的蒙特卡洛算法
1.算法的基本思想
- 连续模型（如在NVT下）
基本思想：
1.首先从一个初始构型开始（原子的位置）
2.进行试探性移动（原子的位移）
3.计算能量变化
4.使用Metropolis算法来接受或拒绝该移动。
2.算法的实现过程
- 算法的实现：
1.初始化：从初始随机构型开始，能量 U n = U o l d U_{n}=U_{old} Un=Uold
2.在所有 N 个原子中随机选取一个原子 i i i，初始位置为 r i , o l d r_{i,old} ri,old
3.试探性移动：通过随机位移将原子 i i i移动到新位置 r i , n e w r_{i,new} ri,new
4.计算能量变化 Δ U = U n e w − U o l d \Delta U = U_{new} - U_{old} ΔU=Unew−Uold
5.根据Metropolis 准则接受或拒绝移动
6.如果试探性移动被接受，保留位移的原子作为新系统 n + 1 n+1 n+1
如果移动被拒绝，新系统n+1保持原来构型
7.计算统计平均值
8.返回至步骤 2 - 我们如何选择合适的位移呢？
蒙特卡洛的位移：
1.首先，位移的方向必须是随机且无偏的；
2.其次，必须选择位移的长度，使体系尽可能高效地在模型空间中移动。
3.最后，由于原子不是在晶格上移动而是连续在空间中运动，所以位移的长度必须是随机的。
这篇关于使用Metropolis蒙特卡洛方法的原子模拟的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
