这里会显示出您选择的修订版和当前版本之间的差别。
两侧同时换到之前的修订记录前一修订版后一修订版 | 前一修订版 | ||
atk:使用metadynamics动力学方法研究cu_111_中cu空位的扩散_plumed [2019/12/28 22:21] – [在 Cu(111) 上创建一个空位] xie.congwei | atk:使用metadynamics动力学方法研究cu_111_中cu空位的扩散_plumed [2020/02/07 15:50] (当前版本) – [分析结果] fermi | ||
---|---|---|---|
行 1: | 行 1: | ||
====== 使用Metadynamics动力学方法研究 Cu(111) 中 Cu 空位的扩散(PLUMED) ====== | ====== 使用Metadynamics动力学方法研究 Cu(111) 中 Cu 空位的扩散(PLUMED) ====== | ||
- | **版本:**P-219.03 | + | **版本:**P-2019.03 |
- | 在本教程中,您将学习如何使用 '' | + | 本教程介绍如何使用 '' |
===== 介绍 ===== | ===== 介绍 ===== | ||
- | Metadynamics <color # | + | Metadynamics <color # |
- | 在 Metadynamics 中,附加的偏势能会在模拟期间定期添加到系统中。这就使系统能够避免较深的自由能最小值并克服高能量位垒,从而有效地探索整个自由能表面。在材料科学领域,该方法已用于研究晶体多态性、固液界面以及固体和表面的化学反应。在这里,您将应用它研究 Cu(111) 面上 Cu 空位的扩散。 | + | 在 Metadynamics |
- | **QuantumATK** 通过 PLUMED 插件实现 Metadynamics。在本教程中,我们将会简要地介绍该方法,解释怎样为 Cu(111) 上的空位扩散设置 Metadynamics 脚本,以及如何分析生成的自由能表面。 | + | **QuantumATK** 通过 PLUMED 插件实现 Metadynamics。在本教程中,我们先简要地介绍该方法,然后演示如何为 Cu(111) 上的空位扩散设置 Metadynamics 脚本,以及如何分析生成的自由能表面。 |
关于 Metadynamics 方法的更多信息,可参见参考文献 <color # | 关于 Metadynamics 方法的更多信息,可参见参考文献 <color # | ||
===== 理论背景 ===== | ===== 理论背景 ===== | ||
- | 在 Metadynamics 中,在几个选定的自由度相空间中构造了外部偏势能,$S$,通常称为变量集(CVs)< | + | 在 Metadynamics 中,在特定的自由度相空间中构造外加的歧视势能,$S$,这些特定的自由度通常称为综合变量(CVs)< |
$$S(R) = (S_1 (R), ..., S_d (R)).$$ | $$S(R) = (S_1 (R), ..., S_d (R)).$$ | ||
- | 势能是在时间 t 内,高斯函数沿 CV 空间中轨迹的积分。在文献 <color # | + | 势能是在时间 t 内,高斯函数沿 CV 空间中轨迹的积分。在文献 <color # |
$$V(S,t) = \int_{0}^{t} dt^{'} \frac{W}{\tau} \exp\left( -\sum_{i=1}^{d} \frac{(S_i (R)-S_i({R}(t^{' | $$V(S,t) = \int_{0}^{t} dt^{'} \frac{W}{\tau} \exp\left( -\sum_{i=1}^{d} \frac{(S_i (R)-S_i({R}(t^{' | ||
- | 这里,$\tau$ 是高斯函数积分的步幅,$\sigma_i$ 是在第 $i$ 个 CV中高斯函数的宽度,$W$ 是高斯函数的高度。偏势能的作用是使系统远离局部最小值,并使其能够访问相空间的新区域。 | + | 这里,$\tau$ 是高斯函数积分的步幅,$\sigma_i$ 是在第 $i$ 个 CV中高斯函数的宽度,$W$ 是高斯函数的高度。歧视势能的作用是使系统远离当前的局部极小值,并使其能够到达相空间的新区域。 |
===== Cu(111) 中 Cu 空位的 Metadynamics 动力学模拟 ===== | ===== Cu(111) 中 Cu 空位的 Metadynamics 动力学模拟 ===== | ||
行 34: | 行 34: | ||
在 {{: | 在 {{: | ||
- | * 设置弥勒指数为 <111> | + | * 设置 |
* 采用下图所示的表面晶格: | * 采用下图所示的表面晶格: | ||
- | {{ : | + | {{ : |
* 设置平面外晶胞矢量为“Non-periodic and normal (slab)” | * 设置平面外晶胞矢量为“Non-periodic and normal (slab)” | ||
* 设置厚度为 6 层,真空间隙为 10 Å。 | * 设置厚度为 6 层,真空间隙为 10 Å。 | ||
- | {{ : | + | {{ : |
接下来,删除最上面那层坐标为 $(x,y) = (0,0)$ 的原子(下图中红色部分),创造出一个空位: | 接下来,删除最上面那层坐标为 $(x,y) = (0,0)$ 的原子(下图中红色部分),创造出一个空位: | ||
- | {{ : | + | {{ : |
最后,点击 Selection Tools {{: | 最后,点击 Selection Tools {{: | ||
- | {{ : | + | {{ : |
===== 创建 Metadynamics 脚本 ===== | ===== 创建 Metadynamics 脚本 ===== | ||
+ | |||
+ | 现在您可以开始创建 metadynamics 模拟的脚本了。将结构从 **Stash** 发送到 {{: | ||
+ | |||
+ | 1.添加以下模块 | ||
+ | |||
+ | * {{: | ||
+ | * {{: | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | 2.双击 {{: | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | 3.打开 {{: | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | 4.最后,点击 **Add Constraints**,固定 Cu(111) 最下面已被标记的那 4 层。 | ||
+ | |||
+ | |||
+ | {{ : | ||
+ | |||
+ | |||
+ | 点击 {{: | ||
+ | |||
+ | |||
+ | <code python> | ||
+ | # ------------------------------------------------------------- | ||
+ | # Calculator | ||
+ | # ------------------------------------------------------------- | ||
+ | |||
+ | potentialSet = EAM_Cu_2001b() | ||
+ | calculator = TremoloXCalculator(parameters=potentialSet) | ||
+ | calculator.setVerletListsDelta(0.25*Angstrom) | ||
+ | |||
+ | bulk_configuration.setCalculator(calculator) | ||
+ | nlprint(bulk_configuration) | ||
+ | bulk_configuration.update() | ||
+ | nlsave(' | ||
+ | |||
+ | # ------------------------------------------------------------- | ||
+ | # Molecular Dynamics | ||
+ | # ------------------------------------------------------------- | ||
+ | |||
+ | initial_velocity = MaxwellBoltzmannDistribution( | ||
+ | temperature=200.0*Kelvin, | ||
+ | remove_center_of_mass_momentum=True | ||
+ | ) | ||
+ | |||
+ | method = Langevin( | ||
+ | time_step=1*femtoSecond, | ||
+ | reservoir_temperature=200*Kelvin, | ||
+ | friction=0.01*femtoSecond**-1, | ||
+ | initial_velocity=initial_velocity, | ||
+ | heating_rate=0*Kelvin/ | ||
+ | ) | ||
+ | |||
+ | fix_atom_indices_0 = [ 0, 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, 35, 36, 37, 38, | ||
+ | 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, | ||
+ | 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63] | ||
+ | |||
+ | constraints = [FixAtomConstraints(fix_atom_indices_0)] | ||
+ | |||
+ | plumed_command = """ | ||
+ | UNITS LENGTH=A TIME=fs ENERGY=96.48533645956869 | ||
+ | p: POSITION ATOM=81 | ||
+ | METAD ARG=p.x,p.y SIGMA=0.025, | ||
+ | PRINT ARG=p.x, | ||
+ | """ | ||
+ | |||
+ | plumed_hook = PlumedMetadynamics( | ||
+ | configuration=bulk_configuration, | ||
+ | timestep=1.0*fs, | ||
+ | plumed_commands=plumed_command, | ||
+ | logfile=' | ||
+ | ) | ||
+ | md_trajectory = MolecularDynamics( | ||
+ | bulk_configuration, | ||
+ | constraints=constraints, | ||
+ | trajectory_filename=' | ||
+ | steps=10000000, | ||
+ | log_interval=1000, | ||
+ | post_step_hook=plumed_hook, | ||
+ | method=method | ||
+ | ) | ||
+ | |||
+ | bulk_configuration = md_trajectory.lastImage() | ||
+ | nlsave(' | ||
+ | </ | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | 您也可以在此处下载输入文件:[[https:// | ||
+ | |||
+ | 在以上输入中,// | ||
+ | 其中有下面四个命令行: | ||
+ | |||
+ | * // | ||
+ | * // | ||
+ | |||
+ | <WRAP center important 100%> | ||
+ | === 注意 === | ||
+ | 这个数值是根据 [[https:// | ||
+ | </ | ||
+ | |||
+ | |||
+ | 我们将使用在每个轴上// | ||
+ | |||
+ | * // | ||
+ | * // | ||
+ | * // | ||
+ | * // | ||
+ | * // | ||
+ | * // | ||
+ | |||
+ | |||
+ | * // | ||
+ | * // | ||
+ | * // | ||
+ | * // | ||
+ | |||
+ | 如果您想了解更多详尽信息,我们推荐您查阅 [[https:// | ||
+ | |||
+ | 完成后,将其发送到 {{: | ||
+ | |||
+ | |||
+ | |||
+ | |||
==== 分析结果 ==== | ==== 分析结果 ==== | ||
- | ===== 参考 ===== | + | 除了标准的 **QuantumATK** 输出,在作业的最后您将获得 '' |
+ | 为画出自由能 $\mathrm{F}$ 关于综合变量 $\mathrm{CV\ 1}$ 和 $\mathrm{CV\ 2}$ 的函数分布图,首选要用 '' | ||
+ | * '' | ||
+ | * '' | ||
+ | * '' | ||
+ | |||
+ | 可以输入以下命令运行脚本 [[https:// | ||
+ | |||
+ | '' | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | |||
+ | 上图显示了自由能分布关于综合变量 $\mathrm{CV\ 1}$ 和 $\mathrm{CV\ 2}$ 的热图。可以看出,在自由能面上有两个极小值,对应于空位从一个位置到相邻位置的扩散。 | ||
+ | |||
+ | 模拟时间内综合变量的演变可以用以下命令运行脚本 [[https:// | ||
+ | |||
+ | '' | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | 生成图展示了 $\mathrm{CV\ 1}$(红色) 和 $\mathrm{CV\ 2}$(蓝色)关于模拟时间的函数演变。可以看出,$\mathrm{CV\ 2}$ 对应于笛卡尔坐标 $y$,在固定值附近($\mathrm{CV\ 2} = 0 \mathrm{Å}$)振荡,因为两个自由能极小值都出现在相同的 $y$ 值处。相反地,$\mathrm{CV\ 1}$ 对应于笛卡尔坐标 $x$,表明模拟势从左边的极小值($\mathrm{CV\ 1} = -2.5 \mathrm{Å}$)开始。第一个极小值被填充直到接近 $2.2$ ns,然后体系移至第二个极小值($\mathrm{CV\ 1} = 0 \mathrm{Å}$)。第二个极小值也被填充后(大约在 6 ns 内)系统在两个极小值间的振荡概率相同,直到模拟结束。 | ||
+ | |||
+ | 自由能势垒也可以通过运行脚本 [[https:// | ||
+ | |||
+ | '' | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | |||
+ | 结果显示势垒高度为 $0.647$ eV,与 <color # | ||
+ | |||
+ | |||
+ | |||
+ | ===== 参考 ===== | ||
+ | * <color # | ||
+ | * [BBP11] (1, 2, 3) Alessandro Barducci, Massimiliano Bonomi, and Michele Parrinello. Metadynamics. W//IREs: Comput. Mol. Sci//., 1(5): | ||
+ | * <color # | ||
+ | * [LP02] (1, 2) Alessandro Laio and Michele Parrinello. Escaping free-energy minima. //Proc. Natl. Acad. Sci//., 99(20): | ||
+ | * 英文原文:https:// |