用户工具

站点工具

本页面的其他翻译:
  • zh

co:pd_100_结构优化

CO/Pd (100)的结构优化

版本:2015.1

在本教程中,您将学习如何用 Fixed 和 Rigid 约束原子并进行结构优化。学习案例为 CO 分子在 Pd(100) 面上的吸附,所选计算工具为 ATK-DFT 计算引擎和标准 GGA 泛函。

具体地,您将:

  1. 优化块体钯晶体;
  2. 切割块体获得 Pd(100) 面,然后弛豫结构;
  3. 表面吸附 CO 分子,按以下 2 步对体系进行弛豫;
    1. 初始弛豫期间,给原子施加 Rigid 约束;
    2. 最终优化时,对底部的 Pd 原子施加 Fixed 约束。
  4. 最后,弛豫 CO 分子,计算吸附能。

提示

本教程使用特定版本的QuantumATK创建,因此涉及的截图和脚本参数可能与您实际使用的版本略有区别,请在学习时务必注意。

钯块体

创建一个新的项目,命名为 “Pd100_CO”,打开 Builder

  • 点击 Add From Database 将钯块体结构导入到 Stash,然后用发送按钮 传送构形到 Script Generator
  • 在 Script Generator 中,双击添加 New Claculator OptimizeGeometry 模块到脚本,将默认输出文件名称改为 Pd_bulk.nc

  • 双击 New Claculator 模块,在随即打开的窗口中设置参数:
    • k-point sampling:9×9×9;
    • exchange-correlation:GGA-PBE 泛函。
  • 打开 OptimizeGeometry 模块,如下图所示设置参数。特别地:
    • 不勾选 Constrain cell 以确保应力弛豫;
    • 降低 Stress tolerance 至 0.005 eV/Å3
    • 点击 Save trajectory 保存弛豫的轨迹,输入轨迹文件名称为 Pd_bulk.nc

提示

通常,保存轨迹文件是一个很好的做法。假如您的结构优化过程中断,可以轻易地从轨迹的最后一步重启计算。如果您想知道如何重启中断的计算,可以参考教程 Restarting stopped calculations

  • 单击 OK,结构优化窗口关闭。
  • 至此,脚本已经完成。点 按钮将其发送到 Job Manager 。保存出现在窗口的脚本并运行计算,耗时数秒钟。如果有需要,您可以在此 ↓Pd_bulk.py 下载脚本。

块体弛豫结束后,计算的输出结果会在 LabFloor 显示。ID 名字含有 gID000gID002 的块体构形文件分别为最初和最终的结构,弛豫轨迹文件 ID 为 gID001

您可以拖拽这些文件到 Viewer 实现结构可视化。在浏览器中,将鼠标停在晶胞边界可以查看晶格常数,也可以在最新生成的 log 文件中查到。

在接下来的计算中,您所用到的都是都是弛豫后的块体结构 gID002

构造 Pd(100) 面并弛豫

  • 拖拽弛豫后的块体构形到 Builder
  • 打开 Builders Surface (Cleave) 插件。
  • 保留默认的 Miller 指数,单击 Next
  • 保留默认的 1×1 表面晶格,单击 Next
  • Out-of-plane cell vector 下方选择 Non-periodic and normalslab),4 层平面金属,顶部和底端分别有 10 Å 和 5 Å 的真空。单击 Finish 创建表面。

  • 发送板状结构到 Script Generator ,脚本中添加 New Claculator OptimizeGeometry TotalEnergy 模块,修改默认输出文件名称为 Pd100.nc

  • 按照以下设置计算器参数:
    • k-point sampling:9×9×1 k 点;
    • exchange-correlation:GGA-PBE;
    • Poisson solver 选择 FFT2D,平板之下设置 Neumann 边界条件,之上用 Dirichlet 边界条件。

提示

这种沿模拟晶胞 C 方向上的特定边界条件设置很适合板状结构的计算,板上的真空在原则上可以扩展到无限远。

  • 打开 OptimizeGeometry,为力的弛豫做出如下设置:
    • 把轨迹保存到 Pd100_traj.nc 文件;
    • 点击 Add Constraints 打开 Constraints Editor;
    • 用鼠标选中 Pd(100) 面最底部的 2 层,点击 Add tag from selection 中。那两个原子现在被标注为 “Selection 0”,然后为其选择 Fixed 约束。

  • TotalEnergy 模块没有需要编辑的,所以到这步脚本已经完成了。保存为 Pd100.py,用 Job Manager 执行计算,大概需要 1 分钟。如果有需要,您可以在这 ↓Pd100.py 下载脚本。

Pd(100) 面的结构优化完成后,结果会立即出现在 LabFloor。您还可以用 Viewer 把弛豫轨迹可视化。

提示

选中 TotalEnergy,点击 Text Representation 插件可以查看弛豫后 Pd(100) 面的总能量,这个能量值在 log 文件中也有输出。

弛豫 CO/Pd(100) 体系

下一步就是在表面上添加 CO 分子,然后弛豫整个体系。您可以下载现成的 CO/Pd(100) 结构 ↓CO_Pd100.py。您还可以通过教程 Building molecule-surface systems:Benzene on Au(111) 学习怎样在表面上添加分子。

正如在介绍中提到的,您需要两步来弛豫 CO/Pd(100) 构形:

  1. 为了快速了解 CO 被吸附的位置,给 CO 施加 Rigid 约束,面施加 Fixed 约束,并设置一个相对较低的电子密度的网格截断。
  2. 最后一步的弛豫不会约束 CO 分子,只固定 Pd(100) 板的底部。

精确弛豫

将提供的 CO/Pd(100) 构形(↓CO_Pd100.py)移动到 Script Generator,添加 New Claculator OptimizeGeometry 模块。输入默认的输出文件名称为 Pd100_CO_rigid.nc,然后编辑脚本:

  • New Claculator 模块,除了将密度网格截断降至 30 Hartree 外,其他参数使用跟之前弛豫 Pd(100) 时相同的设置。

  • OptimizeGeometry 模块的设置也与弛豫 Pd(100) 时相似,只是约束的选择有略微的不同:
    • 给金属板的所有原子加上 Fixed 约束;
    • CO 分子施加 Rigid 约束(把 C 原子和 O 原子添加到组 “Selection 1”)。

注意

在力优化时,被施加了 Rigid 约束的原子团在移动时像刚体。

  • 保存脚本为 Pd100_CO_rigid.py,用 Job Manager 运行。如果串行执行的话大概耗时 5 分钟。您还可以在此 ↓Pd100_CO_rigid.py 下载脚本。

重要

如果您查看弛豫的轨迹会发现,CO 分子稍微偏移了初始位置,但 C-O 键的长度并没有改变,分子也未旋转。这就是刚性约束的作用。

最终弛豫

在此,您将以上一步得到的结果为初始状态,对 CO/Pd(100) 结构进行最终的结构优化。

  • Pd_CO_rigid.nc 里 ID 为 gID002 的结构发送到 Scripter,默认输出文件名设为 Pd_CO.nc
  • 再次添加 New Claculator OptimizeGeometry TotalEnergy 模块,参数编辑与之前步骤中的相似。只是,这次的 ATD-DFT 计算器的网格截断用 75 Hartree,且只对底部的 2 个 Pd 原子加 Fixed 约束。

  • 保存脚本为 Pd_CO.py,在 Job Manager 或终端执行运算。模拟一共需要约 20 分钟,但如果在更多的 CPU 上并行运行的话,过程会快很多。脚本可以在此下载:↓Pd100_CO.py

提示

在计算运行时,您可以用 Viewer 工具监控结构优化的进程。随着弛豫步数的增加,轨迹对象也在持续更新。

计算完成时,您应该看下完整的弛豫轨迹,可以观察到 CO 分子从在 Pd(100) 面上中空的吸附点附近移动到两个 Pd 表面原子间的桥位。

弛豫CO分子

为了计算得到 Pd(100) 面上 CO 的吸附能,您也需要弛豫单独的 CO 分子。我们可以通过以下步骤获得:

  • Builder 里点击 Add New Configuration 添加新的结构;
  • 选择 Molecular Builder,在打开的面板中点击 Fragments Carbon monoxide;

  • 下一步,选择 New Configuration 的 Stash 里的 H 原子,将其替换为 CO 片段,关闭 Molecular Builder 窗口,将 Stash 里文件名称更改为 CO

提示

想要了解更多关于怎样使用 Molecular Builder 的信息,您可以查看教程: Molecular Builder

现在,您只需要弛豫 CO 原子坐标、计算总能量。把结构发送到 Script Generator,添加 New Claculator OptimizeGeometry TotalEnergy 模块。默认输出文件名称处输入 CO.nc,然后按照以下对模块稍作编辑:

  • New Claculator:选择 GGA-PBE 泛函;
  • OptimizeGeometry:您可以保存轨迹文件或减小 force tolerance,但这不是强制性的。
  • 将脚本保存为 CO.py 并运行,计算速度非常快。脚本也可以在此处下载:↓CO.py

您将会发现 C-O 键长的弛豫结果为 1.143 Å,与实验数据非常吻合。

吸附能

现在您已经准备好在完全的单分子层覆盖下计算 CO/Pd(100) 的吸附能,$\Delta E$ 由总能差得到:

$$\Delta E = E_\mathrm{products} - E_\mathrm{reactants} = E_\mathrm{Pd(100)/CO} - E_\mathrm{Pd(100)} - E_\mathrm{CO}$$

您可以用 Text Representations 插件查看这三个能量值,从而得到 $\Delta E$,或者还可以用脚本实现: ↓adsorption_energy.py

计算吸附能时,采用 PBE 和含有 DZP 基组的 FHI 赝势,结果为 $\Delta E$ = -1.97 eV。但是,这个结果被所谓的基组叠加误差影响,所以您应该继续进行以下部分解决这个问题。

平衡校正

基组叠加误差(BSSE)是由 LCAO 基组的不完全造成的,会对不同子体系间的能量差产生重大影响。正如在教程 (BSSE counterpoise correction)中做出的解释,在 CO/Pd(100) 体系中 Pd 和 CO 基组的叠加会增加由人为因素造成的总能量降低,这是因为 Pd(100) 和 CO 子体系可以互相“借”基函数。结果就是吸附能太大。

中和 BSSE 的标准做法是应用所谓的平衡校正。在 ATK 里,通过使用 CounterpoiseCorrected 实现。

因此,您应该计算 CO/Pd(100) 体系的平衡(CP)校正总能量,以便计算得到的 CP 校正吸附能。

$$\Delta E^\mathrm{CP} = E^\mathrm{CP}_\mathrm{Pd(100)/CO} - E_\mathrm{Pd(100)} - E_\mathrm{CO}$$

采用如下所示的脚本(可在此下载:↓Pd_CO_cp.py)。该脚本读取之前弛豫的 CO/Pd(100) 结构和使用的计算器,并创建新的计算器应用在 CP 校正上。然后弛豫结构,计算总能量。

可以用 Job Manager 或终端运行脚本。它应该只执行 4 个非常小的弛豫步骤,在 10 分钟内完成计算。CP 校正会增大 CO/Pd(100) 体系的能量,导致吸附能 $\Delta E^\mathrm{CP}$ = -1.68 eV。

Pd_CO_cp.py
   1   # -------------------------------------------------------------           
   2   # Bulk Configuration                                                      
   3   # -------------------------------------------------------------           
   4                                                                             
   5   bulk_configuration = nlread('Pd100_CO.nc', BulkConfiguration)[-1]         
   6                                                                             
   7   # Add tags                                                                
   8   bulk_configuration.addTags('CO',          [4, 5])                         
   9   bulk_configuration.addTags('Pd',          [0, 1, 2, 3])                   
  10   bulk_configuration.addTags('Selection 0', [0, 1])                         
  11                                                                             
  12   # -------------------------------------------------------------           
  13   # Calculator                                                              
  14   # -------------------------------------------------------------           
  15                                                                             
  16   # Get the calculator settings from non-BSSE calculation.                  
  17   calculator = bulk_configuration.calculator()                              
  18   basis_set = calculator.basisSet()                                         
  19   exchange_correlation = calculator.exchangeCorrelation()                   
  20   numerical_accuracy_parameters = calculator.numericalAccuracyParameters()  
  21   poisson_solver = calculator.poissonSolver()                               
  22                                                                             
  23   # Create BSSE calculator.                                                 
  24   bsse_calculator = counterpoiseCorrected(LCAOCalculator, ["CO", "Pd"])     
  25   calculator = bsse_calculator(                                             
  26    basis_set=basis_set,                                                     
  27    exchange_correlation=exchange_correlation,                               
  28    numerical_accuracy_parameters=numerical_accuracy_parameters,             
  29    poisson_solver=poisson_solver,                                           
  30    )                                                                        
  31                                                                             
  32   bulk_configuration.setCalculator(calculator)                              
  33   nlprint(bulk_configuration)                                               
  34   bulk_configuration.update()                                               
  35   nlsave('Pd100_CO_cp.nc', bulk_configuration)                              
  36                                                                             
  37   # -------------------------------------------------------------           
  38   # Optimize Geometry                                                       
  39   # -------------------------------------------------------------           
  40   constraints = [0, 1]                                                      
  41                                                                             
  42   bulk_configuration = OptimizeGeometry(                                    
  43    bulk_configuration,                                                      
  44    max_forces=0.05*eV/Ang,                                                  
  45    max_steps=200,                                                           
  46    max_step_length=0.2*Ang,                                                 
  47    constraints=constraints,                                                 
  48    trajectory_filename='Pd100_CO_cp.nc',                                    
  49    disable_stress=True,                                                     
  50    optimizer_method=LBFGS(),                                                
  51    )                                                                        
  52   nlsave('Pd100_CO_cp.nc', bulk_configuration)                              
  53   nlprint(bulk_configuration)                                               
  54                                                                             
  55   # -------------------------------------------------------------           
  56   # Total Energy                                                            
  57   # -------------------------------------------------------------           
  58   total_energy = TotalEnergy(bulk_configuration)                            
  59   nlsave('Pd100_CO_cp.nc', total_energy)                                    
  60   nlprint(total_energy)    ''                                                 

参考

co/pd_100_结构优化.txt · 最后更改: 2018/06/15 10:22 由 fermi

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