表面能的计算

f

使用VASP计算离子电导率,似乎比较麻烦

一个能量较低、结构合理的起始三维构型是进行MD模拟的第一步

分子动力学模拟的基本步骤

  • 确定起始构型
  • 选用适当的力场和模拟软件
  • 构建体系和能量最小化
  • 平衡过程
  • 数据采集过程
  • 结果分析

根据系统的传播方式,分子模拟方法可以分为分子动力学MD蒙特卡洛MC两类方法

  • 分子动力学通过数值积分求解牛顿运动方程,生成系统轨迹,可研究结构、动力学和热力学性质。
  • 蒙特卡洛方法基于概率规则生成新构型,侧重采样平衡态结构和热力学性质,不涉及时间动力学。

力场是一组描述分子中原子间相互作用的势能函数参数的集合,用于计算系统的总能量和原子受力,从而驱动分子动力学(MD)模拟或蒙特卡罗(MC)采样:

  • 将量子力学描述的分子相互作用简化为经典力学模型,降低计算成本。
  • 通过参数化的势能函数,预测分子的结构、动力学和热力学性质(如键长、自由能、扩散系数等)。

固定电荷力场:固定原子电荷,不考虑环境引起的电子云变形(极化),计算效率高

  • 生物分子力场
    • AMBER
      • 广泛用于蛋白质、核酸模拟,参数优化针对溶液环境,如 ff14SB、ff19SB
    • CHARMM
      • 广泛用于蛋白质、核酸模拟,参数优化针对溶液环境,如 ff14SB、ff19SB
    • GROMOS
      • 用于膜蛋白和聚合物,参数简洁,计算速度快
  • 小分子力场
    • GAFF
      • 可参数化多种有机小分子,支持自定义分子
    • UFF
      • 覆盖周期表多数元素,适合无机材料和有机分子

极化力场:显示考虑原子极化率(电子云随电场变化),通过附加自由度(如 Drude 振子、诱导偶极)模拟极化效应,精度更高但计算成本显著增加

  • AMOEBA:基于原子极化率和多极矩,适用于极性体系(如水、蛋白质 - 配体相互作用)。
  • Drude 模型:每个原子附加虚拟振子模拟极化,如 CHARMM-Drude 力场。
  • Polarizable AMBER (PAMBS):在 AMBER 框架中引入极化参数。

粗粒化力场(Coarse-Grained Force Fields):将多个原子简化为单个 “珠子”(如蛋白质的每个残基、脂质的疏水尾),降低计算复杂度,适用于大尺度系统(如细胞膜、聚合物)

  • Martini:广泛用于膜蛋白、病毒颗粒和软物质,4:1 粗粒化(4 个原子→1 个珠子)。
  • MARTINI 3.0:支持多尺度模拟,兼容全原子模型。

使用pymatgen提取数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
'''
分析AIMD结果,计算MSD 和 conductivity
'''
import os
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen.core import Structure
from pymatgen.analysis.diffusion.analyzer import DiffusionAnalyzer
import numpy as np
import pickle
import matplotlib.pyplot as plt

# 这一步是读取 XDATCAR,得到一系列结构信息
traj = Trajectory.from_file('XDATCAR')

# 这一步是实例化 DiffusionAnalyzer 的类
# 并用 from_structures 方法初始化这个类; 900 是温度,2 是POTIM 的值,1是间隔步数
# 间隔步数(step_skip)不太容易理解,但是根据官方教程:
# dt = timesteps * self.time_step * self.step_skip

diff = DiffusionAnalyzer.from_structures(traj,'Li',900,2,1)

# 可以用内置的 plot_msd 方法画出 MSD 图像
# 有些终端不能显示图像,这时候可以调用 export_msdt() 方法,得到数据后再自己作图
ax = diff.export_msdt("msd.dat")
#plt.show()

# 接下来直接得到 离子迁移率, 单位是 mS/cm
C = diff.conductivity

with open('result.dat','w') as f:
f.write('# AIMD result for Li-ion\n')
f.write('temp\tconductivity\n')
f.write('%d\t%.2f\n' %(900,C))

画图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import numpy as np
import matplotlib.pyplot as plt

# 读取数据
data = np.loadtxt('msd.txt', comments='#') # 假设你的文件名为 msd.txt

# 提取时间步长和 MSD 数据
time = data[:, 0] # 第一列是时间步长
msd = data[:, 1] # 第二列是总的 MSD

# 绘制 MSD 曲线
plt.figure(figsize=(10, 6))
plt.plot(time, msd, label='MSD', color='blue', linewidth=2)
plt.xlabel('Time (fs)', fontsize=14)
plt.ylabel('Mean Square Displacement (Ų)', fontsize=14)
plt.title('MSD vs Time', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.show()
  1. System preparation
  2. Minimization/Relaxation
  3. Equilibration
  4. Production

系统准备可以说是模拟中最关键的阶段,但在许多情况下却得到了最少的关注;如果你的系统准备有缺陷,这些缺陷可能会证明是致命的。最坏的可能结果是,准备好的系统并非你想要的(例如,它包含不正确的分子或质子化状态),但在化学上却是有效的,并且可以被你的力场很好地描述,因此没有错误地通过了剩余的步骤——实际上,这是系统准备问题中常见的后果。如果系统在随后的平衡步骤中行为良好,不应假设系统已正确准备;这里应该格外小心。

生产流程如下:

生产流程图

为什么需要有恒温器?

分子动力学模拟用于观察和获取某个研究系统的有趣性质。在很多情况下,为了模拟实验室条件下进行的实验,希望从正则(恒温)配分函数中进行抽样。一般来说,如果系统在模拟过程中必须保持温度不变,将采用某种温控算法。

分子动力学模拟的温度通常使用根据等分配定理定义的动能来测量:

32NkBT=<i=1N12mivi2>\frac{3}{2}Nk_BT = <\sum_{i=1}^{N}\frac{1}{2}m_iv_i^2>

带角标的括号表示温度是一个时间平均值。如果我们使用等分配定理来计算分子动力学模拟单时间点的温度快照,而不是时间平均,这个量称为瞬时温度。瞬时温度不一定等于目标温度;实际上,在正则配分函数中,瞬时温度应该围绕目标温度波动。

温控算法通过改变本质上微观可逆的牛顿运动方程来工作。因此,如果希望计算扩散系数等动力学性质,则最好不使用温控器;相反,在系统达到所需温度后,应关闭温控器。然而,虽然所有温控器都会给出非物理动力学,但有些研究发现它们对特定动力学性质的计算影响很小,并且它们在生产模拟中也常被使用

一些恒温器

  • Gaussian
    • 使用了高斯最小约束原理来确定维持瞬时温度所需的最小扰动力。这种温控器不采样正则分布;相反,它采样的是等速(恒动能)配分函数。但是,等速配分函数采样与正则配分函数相同的配置相空间,因此,使用任一配分函数都可以等效地获得位置依赖(结构)平衡性质。然而,速度依赖(动力学)性质在两个配分函数之间不会等效。这种温控器仅适用于某些高级应用。
  • Simple Velocity Rescaling
    • 实施起来最容易的温控器之一,也是最不物理的温控器之一。这种温控器依赖于重新标粒子的动量,使模拟的瞬时温度恰好与目标温度相匹配,不建议使用
  • Berendsen
    • 既不采样正则分布也不采样等速分布,它的使用可能导致模拟伪象,不建议使用
  • Bussi-Donadio-Parrinello (Canonical Sampling Velocity Rescaling)
    • 正确采样了正则配分函数,也被称为“V-rescale”热浴方法,非常常用的热浴方法
  • Andersen
    • 只能用于采样结构性质,因为动力学性质可能会受到突然碰撞的很大影响,不建议使用
  • Langevin
    • 不了解
  • Noés-Hoover
    • 不了解
压浴类型 原理 特点 适用性 局限性
Simple volume rescaling 每次修改系统体积使瞬时压强等于目标压强 操作简单直接 仅可用于模拟初始阶段,帮助系统初步接近目标压强 不能采样正确配分函数,不能用于生产采样;不平滑地接近目标压强,可能导致系统出现不真实的物理问题
Berendsen 将系统耦合到弱相互作用的压力浴,通过缩放因子定期缩放体积 比简单体积缩放更真实地接近目标压强,产生更真实的压力波动 对于平衡开始阶段可能有用 采样配分函数定义不佳,不能保证是 NPT 或 NPH,不应用于生产采样
Andersen 在运动方程中增加额外自由度,将系统耦合到假想压力浴(类似各向同性活塞作用) 采样正确的配分函数 适用于对配分函数采样要求高的场景 是各向同性的,不能对系统的一部分施加各向异性压力
Parrinello - Rahman 对 Andersen 压缩计的扩展,支持模拟盒大小和形状的各向异性缩放 具有 Andersen 压缩计的性质,增加了各向异性支持,适用于固体模拟(可模拟晶体格子中相变形状变化)等 在固体模拟等领域应用优势明显 并非所有模拟器都已实现该压缩计,限制了用户选择
Martyna - Tuckerman - Tobias - Klein (MTTK) 引入替代运动方程,正确采样较小系统配分函数 通常被视为对 Parrinello - Rahman 的改进,适用于小系统 为小系统的模拟提供了更准确的配分函数采样方法 未提及明显独特的局限性(与其他扩展自由度压缩计类似,受模拟器实现情况等限制)

表面能的计算
https://hydrogen1222.xyz/2025/05/01/科研/VASP/Starry Sky/离子电导率的计算/
作者
Storm Talia
发布于
2025年5月1日
许可协议