这是本文档旧的修订版!
版本:2015.1
在本教程中,您将学习如何用Fixed和Rigid约束原子并进行结构优化。学习案例为CO分子在Pd(100)面上的吸附,所选计算工具为ATK-DFT计算引擎和标准GGA泛函。
具体地,您将:
1、优化块体钯晶体;
2、切割块体获得Pd(100)面,然后弛豫结构;
3、表面吸附CO分子,按以下2步对体系进行弛豫;
(1)初始弛豫期间,给原子施加Rigid约束;
(2)最终优化时,对底部的Pd原子施加Fixed约束。
4、最后,弛豫CO分子,计算吸附能。
创建一个新的项目,命名为“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名字含有gID000和gID002的块体构形文件分别为最初和最终的结构,弛豫轨迹文件ID为gID001。
您可以拖拽这些文件到Viewer 实现结构可视化。在浏览器中,将鼠标停在晶胞边界可以查看晶格常数,也可以在最新生成的log文件中查到。
在接下来的计算中,您所用到的都是都是弛豫后的块体结构gID002。
♥ 打开Builders Surface (Cleave) 插件。
♥ 保留默认的Miller指数,单击Next。
♥ 保留默认的1×1表面晶格,单击Next。
♥ 在Out-of-plane cell vector下方选择Non-periodic and normal(slab),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分子,然后弛豫整个体系。您可以下载现成的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(1 0 0)结构进行最终的结构优化。
♥ 把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。
计算完成时,您应该看下完整的弛豫轨迹,可以观察到CO分子从在Pd(100)面上中空的吸附点附近移动到两个Pd表面原子间的桥位。
为了计算得到Pd(1 0 0)面上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,然后按照以下对模块稍作编辑:
♥ 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。
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) ''