费米维基

用户工具

站点工具

您的足迹:
atk:分子动力学基础

差别

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

到此差别页面的链接

两侧同时换到之前的修订记录前一修订版
后一修订版		前一修订版
atk:分子动力学基础 [2016/09/03 17:21] – [设置计算] dong.dongatk:分子动力学基础 [2019/11/16 22:50] (当前版本) – [介绍] dong.dong
行 2: 行 2:
  
 =====介绍===== =====介绍=====
-Virtual NanoLab 和 ATK 可以使分子动力学模拟变得非常简单：只需要向构型（configuration）添加所需的计算器（calculator）和分子动力学模块（MolecularDynamics block），就可以开始计算了。对于不同类型的模拟，什么样的参数和设置是最合适的，在这里给出了一些指导方针。本文简要介绍 VNL-ATK 中分子动力学（MD）功能，并且一步步的解释了如何正确地设置模拟过程以得到想要的结果。 
  
 所谓分子动力学（MD）模拟，是在预先确定条件下（比如温度，压力，应力，外力等等）模拟原子和分子运动的一种方法。分子动力学模拟可以用来研究纳米尺度下的动力学过程，还可以用来计算相图、扩散系数和各种响应函数等诸多性质。 所谓分子动力学（MD）模拟，是在预先确定条件下（比如温度，压力，应力，外力等等）模拟原子和分子运动的一种方法。分子动力学模拟可以用来研究纳米尺度下的动力学过程，还可以用来计算相图、扩散系数和各种响应函数等诸多性质。
-=====方法论=====+ 
 +QuantumATK（和 NanoLab）可以使分子动力学模拟变得非常简单：只需要向构型（configuration）添加所需的计算器（calculator）和分子动力学模块（MolecularDynamics block），就可以开始计算了。对于不同类型的模拟，什么样的参数和设置是最合适的，在这里给出了一些指导方针。本文简要介绍 QuantumATK 中分子动力学（MD）功能，并且一步步的解释了如何正确地设置模拟过程以得到想要的结果。 
 + 
 +<WRAP center info 100%> 
 +=== 提示 === 
 +**本教程使用特定版本的QuantumATK创建，因此涉及的截图和脚本参数可能与您实际使用的版本略有区别，请在学习时务必注意。** 
 +</WRAP> 
 +===== 方法原理 =====
  
 分子动力学模拟本质上是对给定初始构型中原子运动的牛顿方程求数值解。这通常是通过对有限时间步长进行积分计算得到的。原子之间的相互作用（也就是原子间力）可以通过不同的方法（可以是密度泛函理论（DFT）或经典力场）来进行计算。这些力决定了原子的加速度并且使得原子的位置和速度传递到下一个时间步长。多次重复这个过程会产生一系列的构型快照，这些快照描述了系统在相空间中的运动轨迹，从而可以进一步从中分析提取出想要的性质。 分子动力学模拟本质上是对给定初始构型中原子运动的牛顿方程求数值解。这通常是通过对有限时间步长进行积分计算得到的。原子之间的相互作用（也就是原子间力）可以通过不同的方法（可以是密度泛函理论（DFT）或经典力场）来进行计算。这些力决定了原子的加速度并且使得原子的位置和速度传递到下一个时间步长。多次重复这个过程会产生一系列的构型快照，这些快照描述了系统在相空间中的运动轨迹，从而可以进一步从中分析提取出想要的性质。
行 11: 行 17:
 在设置模拟之前你需要选择你所感兴趣的计算类型。要考虑得因素有：总能量是否守恒（孤立体系中的情况）？模拟体系与热浴（heat bath）之间的耦合时是否要保持温度不变？系统是否受到外加压力或者应力？基于这些考虑我们应该选择一套合适的模拟参数：时间步长大小， 积分步长的数量（模拟持续的时间），积分算法，初始温度，约束条件等等。 在设置模拟之前你需要选择你所感兴趣的计算类型。要考虑得因素有：总能量是否守恒（孤立体系中的情况）？模拟体系与热浴（heat bath）之间的耦合时是否要保持温度不变？系统是否受到外加压力或者应力？基于这些考虑我们应该选择一套合适的模拟参数：时间步长大小， 积分步长的数量（模拟持续的时间），积分算法，初始温度，约束条件等等。
  
-类似于真实的实验，进行分子动力学模拟需要一些实证上的知识和经验。在本例中，你将会了解分子动力学模拟的基本组成部分，以及如何使用 ATK 的分子动力学模块来进行分子动力学模拟。+类似于真实的实验，进行分子动力学模拟需要一些实证上的知识和经验。在本例中，你将会了解分子动力学模拟的基本组成部分，以及如何使用 QuantumATK 的分子动力学模块来进行分子动力学模拟。
  
 将这些技巧应用于一个简单、著名的测试体系——块体硅，你将逐步学习到如何选择时间步长大小，如何控制温度和压力，以及如何通过约束条件来固定某些原子。 将这些技巧应用于一个简单、著名的测试体系——块体硅，你将逐步学习到如何选择时间步长大小，如何控制温度和压力，以及如何通过约束条件来固定某些原子。
行 68: 行 74:
 {{ :atk:mdbasical05.png?600 |}} {{ :atk:mdbasical05.png?600 |}}
  
-总体上，系统的能量保持在一个相对于平均值小于10-5量级上的小涨落内。重要的是，平均值是固定的，没有发现漂移。这些发现表明所选分子动力学（MD）的设置较为恰当。+总体上，系统的能量保持在一个相对于平均值小于 $10^{-5}$ 量级上的小涨落内。重要的是，平均值是固定的，没有发现漂移。这些发现表明所选分子动力学（MD）的设置较为恰当。
  
 我们看看如果选择了一个很大的时间步长大小会发生什么。修改初始温度为 5000 K 和时间步长大小为 5 fs，重复上述模拟。你可以再一次使用MolecularDynamics 模块，也可以将 ''md-nve-1fs.py'' 直接拖拽到 **Editor** 按钮。记得改一下轨迹（trajectory）文件名。保存这个新修改的脚本为 ''md-nve-5000K-5fs.py'' 并运行。 我们看看如果选择了一个很大的时间步长大小会发生什么。修改初始温度为 5000 K 和时间步长大小为 5 fs，重复上述模拟。你可以再一次使用MolecularDynamics 模块，也可以将 ''md-nve-1fs.py'' 直接拖拽到 **Editor** 按钮。记得改一下轨迹（trajectory）文件名。保存这个新修改的脚本为 ''md-nve-5000K-5fs.py'' 并运行。
行 76: 行 82:
 {{ :atk:mdbasical06.png?600 |}} {{ :atk:mdbasical06.png?600 |}}
  
-核心要点是，大的时间步长看上去会提高模拟效率，但实际上会增多运动方程数值积分方案的错误。原子力在一个积分步长内近似守恒的基本假设失效。本质上，所选时间步长应该足够小以解析原子的最高振动频率（也就是说，其应该远小于最小的振动周期）。所以如果你有轻原子（比如氢），那么通常你需要选择比只有重原子(比如金)时更小的时间步长。如果在你的计算中有不同的元素，且温度较高或者原子远离其平衡构型（即粒子受到很大的力）时，也需要选择小的时间步长。如果你不知道使用多大的时间步长，对于多数系统你可以选择 1 fs 作为一个比较安全的开始。通过监测感兴趣的条件下的 NVE 模拟总能是否守恒，你可以对大一些的时间步长值进行评估。+核心要点是，大的时间步长看上去会提高模拟效率，但实际上会增加运动方程数值积分方案的误差。原子力在一个积分步长内近似守恒的基本假设失效。本质上，所选时间步长应该足够小以解析原子的最高振动频率（也就是说，其应该远小于最小的振动周期）。所以如果你有轻原子（比如氢），那么通常你需要选择比只有重原子(比如金)时更小的时间步长。如果在你的计算中有不同的元素，且温度较高或者原子远离其平衡构型（即粒子受到很大的力）时，也需要选择小的时间步长。如果你不知道使用多大的时间步长，对于多数系统你可以选择 1 fs 作为一个比较安全的开始。通过监测感兴趣的条件下的 NVE 模拟总能是否守恒，你可以对大一些的时间步长值进行评估。
  
 另一个需要知道的分子动力学模拟的重要方面是，一个系统需要一定的时间来与所设置的外界温度和压力达到平衡。在本例中，系统温度和总能总的涨落的量级达到稳态用了大约 400 fs。这点很重要，因为任何可观察量的测量需要在系统达到平衡之后才可以实现（除非你要研究非平衡现象）。在本例中你需要等待至少 400 fs。总的来说，分子动力学模拟的平衡时间最大可以直到几百纳秒（比如在生物分子或者聚合物系统），并且依赖于系统中出现的最长的弛豫时间。所以你应该仔细检查，所感兴趣的可观察量是否在你的模拟中达到稳态。 另一个需要知道的分子动力学模拟的重要方面是，一个系统需要一定的时间来与所设置的外界温度和压力达到平衡。在本例中，系统温度和总能总的涨落的量级达到稳态用了大约 400 fs。这点很重要，因为任何可观察量的测量需要在系统达到平衡之后才可以实现（除非你要研究非平衡现象）。在本例中你需要等待至少 400 fs。总的来说，分子动力学模拟的平衡时间最大可以直到几百纳秒（比如在生物分子或者聚合物系统），并且依赖于系统中出现的最长的弛豫时间。所以你应该仔细检查，所感兴趣的可观察量是否在你的模拟中达到稳态。
行 87: 行 93:
 你将以室温下结构的平衡态为开始。用上述同样的方法和同样的模块（NewCalculator 和 Optimize ‣ MolecularDynamics) 来准备脚本。实际上，如果你先前的脚本还存储在 **ScriptGenerator** 里，你可以重新使用它。 你将以室温下结构的平衡态为开始。用上述同样的方法和同样的模块（NewCalculator 和 Optimize ‣ MolecularDynamics) 来准备脚本。实际上，如果你先前的脚本还存储在 **ScriptGenerator** 里，你可以重新使用它。
  
-双击 MolecularDynamics 模块打开设置。由于我们要运行 NVT 系综下的模拟，我们这回需要选择一个不同的 Type。对于 NVT 模拟，ATK 提供了四种选项：+双击 MolecularDynamics 模块打开设置。由于我们要运行 NVT 系综下的模拟，我们这回需要选择一个不同的 Type。对于 NVT 模拟，QuantumATK 提供了四种选项：
   - NVT Berendsen((J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak: Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 81, 3684, 1982))    - NVT Berendsen((J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak: Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 81, 3684, 1982)) 
   - NVT Nose Hoover((W. G. Hoover: Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 31, 1695, 1985))   - NVT Nose Hoover((W. G. Hoover: Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 31, 1695, 1985))
行 140: 行 146:
 打开 **MolecularDynamics** 设置。通过查看 Type 下面的选项，你会找到两个运行 NPT 模拟的选项：**NPT Berendsen** 和 **NPT Melchionna** 算法。 打开 **MolecularDynamics** 设置。通过查看 Type 下面的选项，你会找到两个运行 NPT 模拟的选项：**NPT Berendsen** 和 **NPT Melchionna** 算法。
  
-**NPT Berendsen** 方法的工作原理和 NVT 积分器相同，所以不能精确产生正确的物理系综，尽管偏差通常也很小。而且，ATK 中 Berendsen 的实现只考虑压力标量 $P=1/3(P_{xx}+P_{yy}+P_{zz})$,（其中 $P_{ij}=−σ_{ij}$ 为从应力张量 σ 得到的压力张量的分量），视系统为各项同性。如果模拟像液体或者立方晶体这样的各项同性系统，采用标量压力是一个合适的选择。对于各向异性系统，这个方法无法说明不同方向的不同模量。所以，NPT Berendsen 方法应只能被用于各项同性材料的平衡模拟。+**NPT Berendsen** 方法的工作原理和 NVT 积分器相同，所以不能精确产生正确的物理系综，尽管偏差通常也很小。而且，QuantumATK 中 Berendsen 的实现只考虑压力标量 $P=1/3(P_{xx}+P_{yy}+P_{zz})$,（其中 $P_{ij}=−σ_{ij}$ 为从应力张量 σ 得到的压力张量的分量），视系统为各项同性。如果模拟像液体或者立方晶体这样的各项同性系统，采用标量压力是一个合适的选择。对于各向异性系统，这个方法无法说明不同方向的不同模量。所以，NPT Berendsen 方法应只能被用于各项同性材料的平衡模拟。
  
 为了将压力分量 $P_{ij}$ 与固定的压力分别耦合，你更应该选择使用 NPT Melchionna 方法((S. Melchionna, G. Ciccotti, B.L. Holian: Hoover NPT dynamics for systems varying in shape and size. Mol. Phys. 78, 533, 1993))。平衡定温定压 1 bar 下的硅超胞的合适的设置如下图所示。 为了将压力分量 $P_{ij}$ 与固定的压力分别耦合，你更应该选择使用 NPT Melchionna 方法((S. Melchionna, G. Ciccotti, B.L. Holian: Hoover NPT dynamics for systems varying in shape and size. Mol. Phys. 78, 533, 1993))。平衡定温定压 1 bar 下的硅超胞的合适的设置如下图所示。
行 152: 行 158:
 最后，应力耦合编组框允许你选择恒压器（barostat）与哪个压力/应力张量分量进行耦合，这意味着系统只与静水压力/应力耦合，而忽视剪切应力分量。这些设置对当前的例子十分有效。 最后，应力耦合编组框允许你选择恒压器（barostat）与哪个压力/应力张量分量进行耦合，这意味着系统只与静水压力/应力耦合，而忽视剪切应力分量。这些设置对当前的例子十分有效。
  
-要知道， ATK NPT Melchionna 方法总是将所有压力分量当做互相独立的。并且，不仅是单个外界压力，你也可以为每一个分量设定独立的参考应力值来模拟力学剪切或者蠕变实验。ATK参考手册[[http://www.quantumwise.com/documents/manuals/latest/ReferenceManual/index.html/ref.nptmelchionna.html|Reference Manual]]和例子[[http://quantumwise.com/publications/tutorials/item/576-simulating-a-creep-experiment-of-polycrystalline-copper|Simulating a creep experiment of polycrystalline copper]] 为此功能提供了更详细的解释。+要知道， QuantumATK NPT Melchionna 方法总是将所有压力分量当做互相独立的。并且，不仅是单个外界压力，你也可以为每一个分量设定独立的参考应力值来模拟力学剪切或者蠕变实验。QuantumATK参考手册[[http://www.quantumwise.com/documents/manuals/latest/ReferenceManual/index.html/ref.nptmelchionna.html|Reference Manual]]和例子[[http://quantumwise.com/publications/tutorials/item/576-simulating-a-creep-experiment-of-polycrystalline-copper|Simulating a creep experiment of polycrystalline copper]] 为此功能提供了更详细的解释。
 ====编辑脚本==== ====编辑脚本====
-将脚本送到Editor来增加一些额外的功能。为了可视化体积和压力随时间变化的涨落，将这几行代码附在脚本的末尾：+将脚本送到 Editor 来增加一些额外的功能。为了可视化体积和压力随时间变化的涨落，将这几行代码附在脚本的末尾：
 <code python> <code python>
 # Get the number of snapshots in the trajectory # Get the number of snapshots in the trajectory
行 185: 行 191:
  
 ====模拟结果==== ====模拟结果====
-将脚本送到Job Manager来进行模拟。这将会花费几分钟来完成。使用Movie Tool中的trajectory打开轨迹，你就可以监测在模拟过程中晶胞的尺寸和形状的变化。温度在所选的模拟热浴温度300K周围涨落。+将脚本送到 **Job Manager** 来进行模拟。这将会花费几分钟来完成。使用 **Movie Tool** 中的 trajectory 打开轨迹，你就可以监测在模拟过程中晶胞的尺寸和形状的变化。温度在所选的模拟热浴温度300K周围涨落。
  
-你可以打开Fluctuations_NPT_tau100.png来查看体积的涨落。它应该与下图类似。+你可以打开 ''Fluctuations_NPT_tau100.png'' 来查看体积的涨落。它应该与下图类似。
  
 {{ :atk:MDbasical12.png?500 |}} {{ :atk:MDbasical12.png?500 |}}
  
-上面的曲线显示了晶胞体积显著的振荡。这个振荡是Nose-Hoover类的温度和压力控制下的典型结果。只要他们的量级保持在合理的范围内，并且在模拟过程中不会急剧增加，那么就没有必要担心。振荡的周期可以通过设置恒压器（Barostat）时间尺度来改变。这个值选的越小，振荡频率越高。振荡的振幅随着恒压器（barostat）时间尺度的减小而增加。+上面的曲线显示了晶胞体积显著的振荡。这个振荡是 Nose-Hoover 类的温度和压力控制下的典型结果。只要他们的量级保持在合理的范围内，并且在模拟过程中不会急剧增加，那么就没有必要担心。振荡的周期可以通过设置恒压器（Barostat）时间尺度来改变。这个值选的越小，振荡频率越高。振荡的振幅随着恒压器（barostat）时间尺度的减小而增加。
  
-下面的曲线是模拟过程中的压力，显示了类似的行为。它在目标压力1bar附近振荡，但是振幅较大。即时压力值如此大的涨落其实在分子动力学模拟中很正常，因为压力是基于应力来计算的。而应力对于体积变化非常敏感。此外，压力是一个宏观性质，对于如此小的系统考虑相应的平均值比即时值更合理。+下面的曲线是模拟过程中的压力，显示了类似的行为。它在目标压力 1 bar 附近振荡，但是振幅较大。即时压力值如此大的涨落其实在分子动力学模拟中很正常，因为压力是基于应力来计算的。而应力对于体积变化非常敏感。此外，压力是一个宏观性质，对于如此小的系统考虑相应的平均值比即时值更合理。
  
-NPT模拟一个典型的应用是模拟相变，比如对石英在不同温度下进行一系列模拟并测量平均晶胞体积随温度变化的函数，来研究α-石英到β-石英相变过程，如文献((M. H. Müser, and K. Binder: Molecular dynamics study of the α-β transition in quartz: elastic properties, finite size effects, and hysteresis in the local structure. Phys. Chem. Minerals 28, 746, 2001))所报道。+NPT 模拟一个典型的应用是模拟相变，比如对石英在不同温度下进行一系列模拟并测量平均晶胞体积随温度变化的函数，来研究α-石英到β-石英相变过程，如文献((M. H. Müser, and K. Binder: Molecular dynamics study of the α-β transition in quartz: elastic properties, finite size effects, and hysteresis in the local structure. Phys. Chem. Minerals 28, 746, 2001))所报道。
 =====有约束条件的分子动力学模拟===== =====有约束条件的分子动力学模拟=====
-在进行分子动力学模拟过程中你可以固定选中原子的位置。设置这种约束最简单的方式是打开**Molecular Dynamics**工具然后点击 Add constraints 打开Constraints editor 工具。+在进行分子动力学模拟过程中你可以固定选中原子的位置。设置这种约束最简单的方式是打开 **Molecular Dynamics** 工具然后点击 Add constraints 打开 Constraints editor 工具。
  
 {{ :atk:MDbasical13.png?400 }} {{ :atk:MDbasical13.png?400 }}
  
-在工具栏里你会找到一个表格列出了所有被标签标记的原子组。你可以通过在展示窗口中（用鼠标拖一个矩形或者点住 Ctrl键选择多个原子）选择所需原子并点击Add tag from selection来添加新的原子组。可以在Constraint栏中对每一个原子组指定约束条件。第一种将对应原子组里的所有原子固定在它们的初始位置，而第二种将原子组内的原子整体视为一个刚性体。这意味着这些原子根据它们质心所受的力进行协同的移动。在当前的实行中只可能平移运动，而不能进行刚性体旋转。+在工具栏里你会找到一个表格列出了所有被标签标记的原子组。你可以通过在展示窗口中（用鼠标拖一个矩形或者点住 Ctrl 键选择多个原子）选择所需原子并点击 Add tag from selection 来添加新的原子组。可以在 Constraint 栏中对每一个原子组指定约束条件。第一种将对应原子组里的所有原子固定在它们的初始位置，而第二种将原子组内的原子整体视为一个刚性体。这意味着这些原子根据它们质心所受的力进行协同的移动。在当前的实行中只可能平移运动，而不能进行刚性体旋转。
  
-当使用约束时，你必须注意至少让一个原子保持无约束的。像NVT NoseHoover和NPT Melchionna这样的恒温器（thermostat）需要至少两个自由原子，因为它们还固定了系统的质心。恒温器（thermostat）会自动适应由约束所导致的自由度数的减少。+当使用约束时，你必须注意至少让一个原子保持无约束的。像 NVT NoseHoover 和 NPT Melchionna 这样的恒温器（thermostat）需要至少两个自由原子，因为它们还固定了系统的质心。恒温器（thermostat）会自动适应由约束所导致的自由度数的减少。
  
-然而，你需要知道Movie Tool或者MD Analyzer中一些功能不能与约束共同使用。比如，当分子动力学模拟中有约束时，Movie Tool中显示的温度可能不会反应系统自由原子实际的温度，而会显示一个较低温度。你可以通过在脚本里添加动能（kinetic energy）代码来监测温度，如下例所示：+然而，你需要知道 **Movie Tool** 或者 **MD Analyzer** 中一些功能不能与约束共同使用。比如，当分子动力学模拟中有约束时，Movie Tool 中显示的温度可能不会反应系统自由原子实际的温度，而会显示一个较低温度。你可以通过在脚本里添加动能（kinetic energy）代码来监测温度，如下例所示：
  
 <code python> <code python>
行 218: 行 224:
 </code> </code>
 <WRAP center box 100%> <WRAP center box 100%>
-注意+=== 注意 ===
  
-通常不建议在有约束原子的情况下实行定压力/应力（比如NPT）模拟，因为系统约束部分内产生的应力分布没有从全局应力张量中移除。+通常不建议在有约束原子的情况下实行定压力/应力（比如 NPT）模拟，因为系统约束部分内产生的应力分布没有从全局应力张量中移除。
 </WRAP> </WRAP>
  
 =====器件构型===== =====器件构型=====
-在**DeviceConfiguration**中实现NVE/NVT分子动力学模拟本质上是和在**BulkConfiguration**中一样的方式。拓展到中间区域的电极中的原子自动被约束不动以维持与电极的相容性。虽然技术上可以在**DeviceConfiguration**中实现定压力/应力（NPT）模拟，由于约束的存在，不建议这种计算。对于**DeviceConfiguration**晶胞的优化有多种方法供选择（参见例子[[http://quantumwise.com/support/tutorials/item/878-geometry-optimization-of-interfaces|Geometry optimization of interfaces]]）。+在 **DeviceConfiguration** 中实现 NVE/NVT 分子动力学模拟本质上是和在 **BulkConfiguration** 中一样的方式。拓展到中间区域的电极中的原子自动被约束不动以维持与电极的相容性。虽然技术上可以在 **DeviceConfiguration** 中实现定压力/应力（NPT）模拟，由于约束的存在，不建议这种计算。对于 **DeviceConfiguration** 晶胞的优化有多种方法供选择（参见例子[[http://quantumwise.com/support/tutorials/item/878-geometry-optimization-of-interfaces|Geometry optimization of interfaces]]）。 
 + 
 +===== 参考 =====
  
 +  * 原文：[[http://docs.quantumwise.com/tutorials/md_basics.html]]
   * 本文翻译：王吉章   * 本文翻译：王吉章
atk/分子动力学基础.1472894490.txt.gz · 最后更改: 2016/09/03 17:21 由 dong.dong

页面工具

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