费米维基

用户工具

站点工具

您的足迹:
atk:使用metadynamics动力学方法研究cu_111_中cu空位的扩散_plumed

差别

这里会显示出您选择的修订版和当前版本之间的差别。

到此差别页面的链接

两侧同时换到之前的修订记录前一修订版
后一修订版		前一修订版
atk:使用metadynamics动力学方法研究cu_111_中cu空位的扩散_plumed [2019/12/28 22:09] – [理论背景] xie.congweiatk:使用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
  
-在本教程中，您将学习如何使用 ''PlumedMetadynamics'' 类和 [[https://www.plumed.org/|PLUMED]] 插件在 **QuantumATK** 中设置和运行 Metadynamics 模拟。它将涵盖设置探索与 Cu(111) 面上 Cu 空位扩散相关的自由能分布，如何设置 Metadynamics 计算以及如何分析结果。+本教程介绍如何使用 ''PlumedMetadynamics'' 类和 [[https://www.plumed.org/|PLUMED]] 插件在 **QuantumATK** 中设置和运行 Metadynamics 模拟，介绍如何研究 Cu(111) 表面上 Cu 空位扩散相关的自由能分布，以及如何设置 Metadynamics 计算以及如何分析结果。
 ===== 介绍 ===== ===== 介绍 =====
  
  
-Metadynamics <color #00a2e8>[LP02]</color> 是一种基于分子动力学的强大方法，可以研究多维自由能表面（FESs）。它允许人们通过基于系统微观坐标的多个变量集（CV）的函数重构系统的 FES。+Metadynamics <color #00a2e8>[LP02]</color> 是基于分子动力学的一种很有用的模拟方法，可以研究多维自由能表面（FESs），通过基于体系微观坐标的多个综合变量（Collective Variables，CV）的函数重构体系的 FES。
  
-在 Metadynamics 中，附加的偏势能会在模拟期间定期添加到系统中。这就使系统能够避免较深的自由能最小值并克服高能量位垒，从而有效地探索整个自由能表面。在材料科学领域，该方法已用于研究晶体多态性、固液界面以及固体和表面的化学反应。在这里，您将应用它研究 Cu(111) 面上 Cu 空位的扩散。+在 Metadynamics 模拟过程中，会在模拟期间定期向体系中添加附加的“歧视”势能（bias potential）。这可以使体系能够克服很高势垒，逃离较深的自由能极小值，从而有效地探索整个自由能表面。在材料科学领域，该方法已用于研究晶体多态性、固液界面以及固体和表面的化学反应。本教程将应用这种方法研究 Cu(111) 表面上 Cu 空位的扩散。
  
-**QuantumATK** 通过 PLUMED 插件实现 Metadynamics。在本教程中，我们将会简要地介绍该方法，解释怎样为 Cu(111) 上的空位扩散设置 Metadynamics 脚本，以及如何分析生成的自由能表面。+**QuantumATK** 通过 PLUMED 插件实现 Metadynamics。在本教程中，我们先简要地介绍该方法，然后演示如何为 Cu(111) 上的空位扩散设置 Metadynamics 脚本，以及如何分析生成的自由能表面。
  
 关于 Metadynamics 方法的更多信息，可参见参考文献 <color #00a2e8>[BBP11]</color> 和 <color #00a2e8>[BN15]</color>。 关于 Metadynamics 方法的更多信息，可参见参考文献 <color #00a2e8>[BBP11]</color> 和 <color #00a2e8>[BN15]</color>。
 ===== 理论背景 ===== ===== 理论背景 =====
  
-在 Metadynamics 中，在几个选定的自由度相空间中构造了外部偏势能，$S$，通常称为变量集（CVs）<color #00a2e8>[LP02]</color>。$S$ 是一系列 d 个关于系统微观坐标 $R$ 的函数。+在 Metadynamics 中，在特定的自由度相空间中构造外加的歧视势能，$S$，这些特定的自由度通常称为综合变量（CVs）<color #00a2e8>[LP02]</color>。$S$ 是一系列 d 个体系微观坐标 $R$ 的函数。
  
 $$S(R) = (S_1 (R), ..., S_d (R)).$$ $$S(R) = (S_1 (R), ..., S_d (R)).$$
  
-势能是在时间 t 内，高斯函数沿 CV 空间中轨迹的积分。在文献 <color #00a2e8>[BBP11]</color> 中，势能公式的简化版本为：+势能是在时间 t 内，高斯函数沿 CV 空间中轨迹的积分。在文献 <color #00a2e8>[BBP11]</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^{'} )))^2}{2\sigma_i^2} \right),$$ $$V(S,t) = \int_{0}^{t} dt^{'} \frac{W}{\tau} \exp\left( -\sum_{i=1}^{d} \frac{(S_i (R)-S_i({R}(t^{'} )))^2}{2\sigma_i^2} \right),$$
  
-这里，$\tau$ 是高斯函数积分的步幅，$\sigma_i$ 是在第 $i$ 个 CV中高斯函数的宽度，$W$ 是高斯函数的高度。偏势能的作用是使系统远离局部最小值，并使其能够访问相空间的新区域。+这里，$\tau$ 是高斯函数积分的步幅，$\sigma_i$ 是在第 $i$ 个 CV中高斯函数的宽度，$W$ 是高斯函数的高度。歧视势能的作用是使系统远离当前的局部极小值，并使其能够到达相空间的新区域。
 ===== Cu(111) 中 Cu 空位的 Metadynamics 动力学模拟 ===== ===== Cu(111) 中 Cu 空位的 Metadynamics 动力学模拟 =====
  
行 30: 行 30:
 ==== 在 Cu(111) 上创建一个空位 ==== ==== 在 Cu(111) 上创建一个空位 ====
  
 +想要在 Cu(111) 上构造一个 Cu 空位，首先您需要创建一个  Cu(111) 面。
 +
 +在 {{:atk:builder.png?direct&25|}} **Builder**，点击 Add {{:atk:arrow.png?direct&5|}} From Database，导入“//Copper//”。然后点击 Builders {{:atk:arrow.png?direct&5|}} Surface (Cleave)，利用以下设置构造一个 Cu(111) 面：
 +
 +  * 设置 Miller 指数为 <111>
 +  * 采用下图所示的表面晶格：
 +
 +{{ :atk:surface_lattice_424-20191228.png?500 |}}
 +
 +  * 设置平面外晶胞矢量为“Non-periodic and normal (slab)”
 +  * 设置厚度为 6 层，真空间隙为 10 Å。
 +
 +{{ :atk:surface_slab-20191228.png?500 |}}
 +
 +接下来，删除最上面那层坐标为 $(x,y) = (0,0)$ 的原子（下图中红色部分），创造出一个空位：
 +
 +{{ :atk:vacancy_site-20191228.png?500 |}}
 +
 +
 +最后，点击 Selection Tools {{:atk:arrow.png?direct&5|}} Tags，标记 Cu(111) 最下面的 4 层。此标记将用于在 metadynamics 模拟过程中对它们的固定。
 +
 +{{ :atk:tag_selection0-20191228.png?500 |}}
  
 ===== 创建 Metadynamics 脚本 ===== ===== 创建 Metadynamics 脚本 =====
 +
 +现在您可以开始创建 metadynamics 模拟的脚本了。将结构从 **Stash** 发送到 {{:atk:script_generator.png?direct&25|}} **Script Generator** ，按照如下设置脚本：
 +
 +1.添加以下模块
 +
 +  * {{:atk:calculator.png?direct&25|}} Calculators {{:atk:arrow.png?direct&5|}} ForceFieldCalculator
 +  * {{:atk:optimization.png?direct&25|}} Optimization {{:atk:arrow.png?direct&5|}} MolecularDynamics
 +
 +{{ :atk:script2-20191228.png?500 |}}
 +
 +2.双击 {{:atk:calculator.png?direct&25|}} **ForceFieldCalculator** 模块，使用“//EAM_Cu_2001b//”参数组。
 +
 +{{ :atk:new_calculator2-20191228.png?500 |}}
 +
 +3.打开 {{:atk:optimization.png?direct&25|}} **MolecularDynamics** 模块，如下图设置模拟参数：
 +
 +{{ :atk:script_md-20191228.png?500 |}}
 +
 +4.最后，点击 **Add Constraints**，固定 Cu(111) 最下面已被标记的那 4 层。
 +
 +
 +{{ :atk:md_constraints_editor-20191228.png?500 |}}
 +
 +
 +点击 {{:atk:script_generator.png?direct&25|}} **Scripter** 里的 {{:atk:sendto.png?direct&20|}} 按钮，将脚本发送到 {{:atk:editor.png?direct&25|}} **Editor**。然后添加下文所示高亮部分，包含了执行 metadynamics 计算的 ''[[https://docs.quantumatk.com/manual/Types/PlumedMetadynamics/PlumedMetadynamics.html#NL.Dynamics.MolecularDynamics.PlumedMetadynamics.PlumedMetadynamics|PlumedMetadynamics]]'' 类。
 +
 +
 +<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('vacancy-hopping.hdf5', bulk_configuration)
 +
 +# -------------------------------------------------------------
 +# 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/picoSecond,
 +)
 +
 +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,0.025 HEIGHT=0.05,0.05 PACE=1000 LABEL=restraint
 +PRINT ARG=p.x,p.y,restraint.bias STRIDE=10 FILE=COLVAR
 +"""
 +
 +plumed_hook = PlumedMetadynamics(
 +    configuration=bulk_configuration,
 +    timestep=1.0*fs,
 +    plumed_commands=plumed_command,
 +    logfile='plumed.log',
 +)
 +md_trajectory = MolecularDynamics(
 +    bulk_configuration,
 +    constraints=constraints,
 +    trajectory_filename='vacancy-hopping-trajectory.hdf5',
 +    steps=10000000,
 +    log_interval=1000,
 +    post_step_hook=plumed_hook,
 +    method=method
 +)
 +
 +bulk_configuration = md_trajectory.lastImage()
 +nlsave('vacancy-hopping.hdf5', md_trajectory)
 +</code>
 +
 +{{ :atk:代码高亮-20191228.png |}}
 +
 +您也可以在此处下载输入文件：[[https://docs.quantumatk.com/_downloads/vacancy-hopping.py|↓ vacancy-hopping.py]]。
 +
 +在以上输入中，//plumed_command// 部分包括被直接传递到 PLUMED <color #00a2e8>[BBP11]</color> 的输入。
 +其中有下面四个命令行：
 +
 +  * //UNITS//：这行控制 PLUMED 中采用的单位。如上编辑此行，将会用到 //Å//、//fs// 和 //eV// 这些单位。
 +  * //p//：这行是用来定义约束的原子，在本例中是序号为 $81$ 的原子。
 +
 +<WRAP center important 100%>
 +=== 注意 ===
 +这个数值是根据 [[https://www.plumed.org/|PLUMED]] 中使用的序号设置的，在 **QuantumATK** 中原子的序号从 $1$ 开始，而不是 $0$。因此这个值应被设置为“//序号//+1 号原子”，这里的原子//序号//是指被约束原子的序号。
 +</WRAP>
 +
 +
 +我们将使用在每个轴上//宽度// $0.025$ 、//高度// $0.05$ 的高斯函数查看 $x$ 轴和 $y$ 轴上的势能面。
 +
 +  * //METAD//：这行控制模拟期间 metadynamics 的类型和采用的高斯函数。
 +      * //ARG//：这个参数用于设置变量集。在本例中，我们将使用的变量集为原子 $p$ 沿笛卡尔坐标 $x\ (p.x)$ 和 $y\ (p.y)$ 方向。
 +      * //SIGMA//：这个参数用于设置每个变量集所用高斯函数 $\sigma$ 的宽度。本例中，我们为两个变量集都设置为 $\sigma = 0.025$。
 +      * //HEIGHT//：这个参数用于设置每个变量集所用高斯函数的高度。本例中，我们为两个变量集都设置为 $0.05$ eV。
 +      * //PACE//：这个参数用于设置添加新高斯函数的间隔。本例我们每 $1000$ MD 步添加一个新的高斯函数。
 +      * //LABEL//：这个参数用于设置变量集的名称。
 +
 +
 +  * //PRINT//：这行控制 <color #00a2e8>PLUMED</color> 输出的格式。
 +      * //ARG//：这个参数控制输出的变量，本例中为 $p.x$、$p.y$ 和偏差。
 +      * //STRIDE//：表明输出间隔，本例为 10 MD 步。
 +      * //FILE//：此行指明输出文件名称。
 +
 +如果您想了解更多详尽信息，我们推荐您查阅 [[https://www.plumed.org/doc-v2.3/user-doc/html/index.html|PLUMED documentation]]。这些命令都被转至 PlumedMetadynamics 类，后面的使用 **hook 函数**传至 ''[[https://docs.quantumatk.com/manual/Types/PlumedMetadynamics/PlumedMetadynamics.html#NL.Dynamics.MolecularDynamics.PlumedMetadynamics.PlumedMetadynamics|MolecularDynamics]]'' 类。
 +
 +完成后，将其发送到 {{:atk:job_manager.png?direct&25|}} **Job Manager**，运行模拟，在 8 个 CPU 的情况下需要 8 个小时。
 +
 +
 +
 +
  
  
 ==== 分析结果 ==== ==== 分析结果 ====
  
-===== 参考 =====+除了标准的 **QuantumATK** 输出，在作业的最后您将获得 ''COLVAR'' 和 ''HILLS'' 文件。这两个 <color #00a2e8>PLUMED</color> 的输出文件将用于分析 metadynamics 模拟。
  
 +为画出自由能 $\mathrm{F}$ 关于综合变量 $\mathrm{CV\ 1}$ 和 $\mathrm{CV\ 2}$ 的函数分布图，首选要用 ''atkpython extract_F.py'' 命令运行脚本 [[https://docs.quantumatk.com/_downloads/extract_F.py|↓ extract_F.py]]。将会产生这三个输出文件：
  
 +  * ''F_vs_cv1.dat'' 包含 $\mathrm{F}$ vs. $\mathrm{CV\ 1}$ 的数据。 
 +  * ''F_vs_cv2.dat'' 包含 $\mathrm{F}$ vs. $\mathrm{CV\ 2}$ 的数据。
 +  * ''F_vs_cv1_cv2.dat'' 包含 $\mathrm{F}$ vs. $\mathrm{CV\ 1}$ 和 $\mathrm{F}$ vs. $\mathrm{CV\ 2}$ 的数据。
 +
 +可以输入以下命令运行脚本 [[https://docs.quantumatk.com/_downloads/plot_F_vs_cv1_cv2.py|↓ plot_F_vs_cv1_cv2.py]] 即可得到自由能分布 $\mathrm{F}$ 的示意图：
 +
 +''atkpython F_vs_cv1_cv2.py''
 +
 +{{ :atk:f_vs_cv1_cv2-20191228.png?500 |}}
 +
 +
 +上图显示了自由能分布关于综合变量 $\mathrm{CV\ 1}$ 和 $\mathrm{CV\ 2}$ 的热图。可以看出，在自由能面上有两个极小值，对应于空位从一个位置到相邻位置的扩散。
 +
 +模拟时间内综合变量的演变可以用以下命令运行脚本 [[https://docs.quantumatk.com/_downloads/cv1_cv2_vs_time.py|↓ cv1_cv2_vs_time.py]] 获得：
 +
 +''atkpython cv1_cv2_vs_time.py''
 +
 +{{ :atk:cv1_cv2_vs_time-20191228.png?600 |}}
 +
 +生成图展示了 $\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://docs.quantumatk.com/_downloads/barrier.py|↓ barrier.py]] 分析：
 +
 +''atkpython barrier.py''
 +
 +{{ :atk:barrier-20191228.png?600 |}}
 +
 +
 +结果显示势垒高度为 $0.647$ eV，与 <color #00a2e8>[KNB17]</color> 中结果一致。
 +
 +
 +
 +===== 参考 =====
  
  
 +  * <color #00a2e8>[BN15]</color> Kristof M. Bal and Erik C. Neyts. Merging metadynamics into hyperdynamics: Accelerated molecular simulations reaching time scales from microseconds to seconds. //Journal of Chemical Theory and Computation//, 11(10):4545–4554, 2015.
 +  * [BBP11] (1, 2, 3) Alessandro Barducci, Massimiliano Bonomi, and Michele Parrinello. Metadynamics. W//IREs: Comput. Mol. Sci//., 1(5):826–843, 2011. [[http://dx.doi.org/10.1002/wcms.31|doi:10.1002/wcms.31]].
 +  * <color #00a2e8>[KNB17]</color> Suresh Kondati Natarajan and Jörg Behler. Self-diffusion of surface defects at copper-water interfaces. //The Journal of Physical Chemistry C//, 121(8):4368–4383, 2017.
 +  * [LP02] (1, 2) Alessandro Laio and Michele Parrinello. Escaping free-energy minima. //Proc. Natl. Acad. Sci//., 99(20):12562–12566, 2002.
 +  * 英文原文：https://docs.quantumatk.com/tutorials/metadynamics_with_plumed/metadynamics_with_plumed.html
atk/使用metadynamics动力学方法研究cu_111_中cu空位的扩散_plumed.1577542190.txt.gz · 最后更改: 2019/12/28 22:09 由 xie.congwei

页面工具

© 2014-2022 费米科技（京ICP备14023855号）