24小时热门版块排行榜&&&&
当前主题已经存档。
【有奖交流】积极回复本帖子，参与交流，就有机会分得作者 阿勇 的 4 个金币
(小有名气)
在线: 16.5小时
虫号: 33779
注册: 性别: GG专业: 凝聚态物性 II ：电子结构
【求助】VASP中 原子能级问题
想用VASP计算得到H，O，Zn各个原子轨道能级，想用他们的轨道能级来分析他们的成键.
& &首先，我分别在20*20*20 的大的supecell里 放入了单个原子，k点选G点，计算得到他们的轨道能量分别是:
&&H:& && && &&&k-point& &1 :& && & 0.0000& & 0.0000& & 0.0000
& && && && && & band No.&&band energies& &&&occupation
& && && && && && &1& && &-6.3540& && &1.00000& & (H: 1s)
& && && && && && &2& && &-0.0494& && &0.00000
& && && && && && &3& && & 0.2821& && &0.00000
& && && && && && &4& && & 0.2918& && &0.00000
O:& && && && &k-point& &1 :& && & 0.0000& & 0.0000& & 0.0000
& && && && && & band No.&&band energies& &&&occupation
& && && && && & 1& &&&-23.7275& && &2.00000& && & (O:2s)
& && && && && & 2& && &-9.1859& && &1.33333& && & (O:2p)
& && && && && & 3& && &-9.1859& && &1.33333& && & (O:2p)
& && && && && & 4& && &-9.1859& && &1.33333& && &&&O:2p)
& && && && && & 5& && &-0.1450& && &0.00000
& && && && && & 6& && & 0.1190& && &0.00000
& && && && && & 7& && & 0.3438& && &0.00000
& && && && && & 8& && & 0.3443& && &0.00000
& && && && && & k-point& &1 :& && & 0.0000& & 0.0000& & 0.0000
& && && && && &band No.&&band energies& &&&occupation
& && && && && &1& &&&-10.4277& && &2.00000& &&& (Zn:3d)
& && && && && &2& &&&-10.4277& && &2.00000& && &(Zn:3d)
& && && && && &3& &&&-10.4277& && &2.00000& &&& (Zn:3d)& && && && &&&
& && && && && &4& &&&-10.4277& && &2.00000& && &(Zn:3d)
& && && && && &5& &&&-10.4277& && &2.00000& &&& (Zn:3d)
& && && && && &6& && &-6.1993& && &2.00000& && & (Zn:4s)
& && && && && & 7& && &-1.2219& && &0.00000
& && && && && & 8& && &-1.2219& && &0.00000
& && && && && &9& && &-1.2219& && &0.00000
& & 我想问的是把H，O,Zn 的各个能级放在一起比较可不可行，也就是说，由VASP计算得到的各个原子对应的能量他们的参考零点是否相同？
我还把O,Zn原子放在同一个supercell 里计算，O Zn之间的距离大概有17A，认为这样得到的轨道能级能反映出他们各自的原子轨道能级，而其他们轨道能量的参考标准应该相同。但计算完后的结果有些不可思议：
& & O, Zn& &: k-point& &1 :& && & 0.0000& & 0.0000& & 0.0000
&&band No.&&band energies& &&&occupation
& && &1& &&&-21.6514& && &2.00000& & (O:2s)
& && &2& &&&-11.7000& && &2.00000& & (Zn:3d)
& && &3& &&&-11.7000& && &2.00000& & (Zn:3d)
& && &4& &&&-11.6998& && &2.00000& & (Zn:3d)& &&&
& && &5& &&&-11.6998& && &2.00000& & (Zn:3d)& && &
& && &6& &&&-11.6998& && &2.00000& & (Zn:3d)& && &
& &&&7& && &-7.2443& && &1.829140&&&&(Zn:4s)
& && &8& && &-7.1836& && &1.39029& &&& (O:2p)
& && &9& && &-7.1836& && &1.39029& &&& (O:2p)
& &&&10& && &-7.1836& && &1.39029& &
(O:2p)& &&&
& &&&11& && &-2.0005& && &0.00000
& &&&12& && &-2.0005& && &0.00000
& &&&13& && &-2.0005& && &0.00000
& && &&&....
结果发现 Zn 4s 与O 2p之间有电荷转移, 两个原子相隔这么远还能有电荷转移？
后来我把胞扩大到30*30*30 后发现还是只不过转移量小了一些而已
(正式写手)
在线: 231.4小时
虫号: 456224
注册: 专业: 凝聚态物性 II ：电子结构
★ ★ ★ 阿勇(金币+2,VIP+0):谢谢 1-11 13:14aylayl08(金币+1,VIP+0):谢谢解答 1-11 19:24
第一个问题，分别看一下三个原子所在真空区的静电势，得到真空能级，并以此为尺度。
第二个问题，不知道，仔细查一下坐标文件。
[ Last edited by westmonster on
at 11:40 ]
蝉，于其脱壳展翅，则蛰居地下，似无声无息，实则以备有所为。
(小有名气)
在线: 149.1小时
虫号: 934744
注册: 性别: GG专业: 凝聚态物性 II ：电子结构
★ ★ ★ ★ ★ ★ ★ 阿勇(金币+2,VIP+0):先谢谢了 1-11 13:15aylayl08(金币+5,VIP+0):谢谢详细解释 1-11 19:22
我认为你的计算是有问题的！
首先，VASP是采用平面波基组而非LCAO，所以根本不要去考虑什么mulliken布居分布，根本不存在。电荷转移只能表征在空间上，而不能说是在哪个原子上面。空间上的转移可以通过计算电荷密度差来研究，原子布居分部目前只能用bader势来看。
然后，你看到的所谓电荷转移，其实是出现了轨道分数占据态。能带中不应该出现分数占据，出现分数坐标说明你的计算中使用了smearing方法来加速SCF收敛，比如说ISMEAR=2，这个时候虽然收敛的很好很速度，但是得到的结果根本不是能量最低。结构已然不对怎么能算对电子结构呢？要算能带结构，ISMEAR要设置为-5才行！
最后，你的原子距离调整到很大，体系的总轨道就是两原子轨道的简单叠加，可以分别算一个各个原子的能级，分开来看比较好。
[ Last edited by happyjwx on
at 12:30 ]
夜半追忆往昔
(小有名气)
在线: 16.5小时
虫号: 33779
注册: 性别: GG专业: 凝聚态物性 II ：电子结构
引用回帖:Originally posted by westmonster at
第一个问题，分别看一下三个原子所在真空区的静电势，得到真空能级，并以此为尺度。
第二个问题，不知道，仔细查一下坐标文件。
[ Last edited by westmonster on
at 11:40 ] 不好意思，请问怎么看“三个原子所在真空区的静电势，得到真空能级，并以此为尺度”啊？vasp 输出文件中哪个有这个静电势啊？
(小有名气)
在线: 149.1小时
虫号: 934744
注册: 性别: GG专业: 凝聚态物性 II ：电子结构
引用回帖:Originally posted by 阿勇 at
不好意思，请问怎么看“三个原子所在真空区的静电势，得到真空能级，并以此为尺度”啊？vasp 输出文件中哪个有这个静电势啊？ 静电势？LOCPOT文件里面有
但是求静点势有什么意义？
我觉得楼上的意思应该是取一个相同大小的晶胞，不放原子，直接计算真空能，再从总能里面减去。不过好像没什么意义嘛。
夜半追忆往昔
(小有名气)
在线: 16.5小时
虫号: 33779
注册: 性别: GG专业: 凝聚态物性 II ：电子结构
引用回帖:Originally posted by happyjwx at
我认为你的计算是有问题的！
首先，VASP是采用平面波基组而非LCAO，所以根本不要去考虑什么mulliken布居分布，根本不存在。电荷转移只能表征在空间上，而不能说是在哪个原子上面。空间上的转移可以通过计算电荷密 ... 我用的是ISMEAR=0算的，我其实也不是想算能带，只是想看看各个原子轨道的相对位置，以为在同一个胞里，他们的参考标准会一样，就不用做真空能级的修正了（其实是我不知道怎么做）
(小有名气)
在线: 149.1小时
虫号: 934744
注册: 性别: GG专业: 凝聚态物性 II ：电子结构
★ ★ ★ ★ aylayl08(金币+4,VIP+0):谢谢指点 1-11 19:23
就算是不算能带的话，你的结构本身可能也有问题。
如果要比较好的收敛，那么ISMEAR=-3
这样的话可以多设置一些SMEAR能量值，从大到小，一直到收敛为止
比如说0.05 0.04 0.03 0.02 0.01 0.00
这样从0.05到0.00达到真正的收敛
可能有些麻烦，如果感兴趣的话你可以试一下
同一个胞里面一般是认为一样的，这一点一般没有人去怀疑吧
如果这都不一样的话，就没什么比较的办法了
夜半追忆往昔
(小有名气)
在线: 16.5小时
虫号: 33779
注册: 性别: GG专业: 凝聚态物性 II ：电子结构
引用回帖:Originally posted by happyjwx at
就算是不算能带的话，你的结构本身可能也有问题。
如果要比较好的收敛，那么ISMEAR=-3
这样的话可以多设置一些SMEAR能量值，从大到小，一直到收敛为止
比如说0.05 0.04 0.03 0.02 0.01 0.00
这样从0.05到0.00 ... 好的，我先试试,谢谢了...
(小有名气)
在线: 16.5小时
虫号: 33779
注册: 性别: GG专业: 凝聚态物性 II ：电子结构
引用回帖:Originally posted by happyjwx at
就算是不算能带的话，你的结构本身可能也有问题。
如果要比较好的收敛，那么ISMEAR=-3
这样的话可以多设置一些SMEAR能量值，从大到小，一直到收敛为止
比如说0.05 0.04 0.03 0.02 0.01 0.00
这样从0.05到0.00 ... ISMEAR=-3 怎么设置SIGMA啊？ 在VASP手册上有：
“SMEARINGS= ismear1 sigma1 ismear2 sigma2 ...”
这是什么意思？
我用如下设置：
ISMEAR = -3
SIGMA = 0.05 0.04 0.03 0.02 0.01
出现错误如下：
internal error: XML_INCAR called with a vector, use XML_INCAR_V instead SIGMA
(小有名气)
在线: 149.1小时
虫号: 934744
注册: 性别: GG专业: 凝聚态物性 II ：电子结构
★ ★ ★ 阿勇(金币+2,VIP+0): 1-11 22:22zxzj05(金币+1,VIP+0):给个红包，谢谢回帖交流 1-12 21:26
引用回帖:Originally posted by 阿勇 at
ISMEAR=-3 怎么设置SIGMA啊？ 在VASP手册上有：
“SMEARINGS= ismear1 sigma1 ismear2 sigma2 ...”
这是什么意思？
我用如下设置：
ISMEAR = -3
SIGMA = 0.05 0.04 0.03 0.02 0.01
出现错误如下：
int ... 按照说明书上写啊
SMEARINGS=0 0.05 0 0.04 0 0.03 0 0.02 0 0.01 0 0.00
这样写。一个ISMEAR一个SIGMA值
夜半追忆往昔
相关版块跳转
第一性原理
我要订阅楼主
的主题更新
小木虫,学术科研互动社区,为中国学术科研免费提供动力
违规贴举报删除请联系客服电话： 邮箱：（全天候） 或者 QQ:
广告投放与宣传请联系 李想 QQ：
QQ：&&邮箱：
Copyright & 2001-, All Rights Reserved. 小木虫 版权所有
浏览器进程
打开微信扫一扫
随时随地聊科研当前位置： >>
第5章-第一性原理计算-2-CASTEP
第五章第一性原理计算 First-principles Calculations讲授：胡朝浩 材料科学与工程学院 桂林电子科技大学第一性原理计算学习和使用 - CASTEP桂林电子科技大学http://www.<b
r />关于CASTAPCASTAP是特别为固体材料学而设计的一个现代的量子力学基本 程序，其使用了密度泛函(DFT）平面波赝势方法，进行第一原理量 子力学计算，以探索如半导体，陶瓷，金属，矿物和沸石等材料的 晶体和表面性质。 典型的应用包括表面化学，键结构，态密度和光学性质等研究， CASTAP也可用于研究体系的电贺密度和波函数的3D形式。此外， CASTAP可用于有效研究点缺陷（空位，间隙和置换杂质）和扩展 缺陷（如晶界和位错）的性质。 Material Studio使用组件对话框中的CASTAP选项允许准备，启 动，分析和监测CASTAP服役工作。 计算：允许选择计算选项（如基集，交换关联势和收敛判据），作 业控制和文档控制。 分析：允许处理和演示CASTAP计算结果。这一工具提供加速整体 直观化以及键结构图，态密度图形和光学性质图形。CASTAP的任务CASTAP计算是要进行的三个任务中的一个，即单个点的能量 计算，几何优化或分子动力学。可提供这些计算中的每一个以便产 生特定的物理性能。性质为一种附加的任务，允许重新开始已完成 的计算以便产生最初没有提出的额外性能。 在CASTAP计算中有很多运行步骤，可分为如下几组： ? 结构定义：必须规定包含所感兴趣结构的周期性的3D模型文 件，有大量方法规定一种结构：可使用构建晶体（Build Crystal)或 构建真空板(Build Vacuum Stab)来构建，也可从已经存在的的结构 文档中引入，还可修正已存在的结构。注意： CASTAP仅能在3D周期模型文件基础上进行计算，必须构建超单胞，以便研 究分子体系。 提示： CASTAP计算所需时间随原子数平方的增加而增加。因此，建议是用最小的初 晶胞来描述体系，可使用Build\Symmetry\Primitive Cell菜单选项来转换成初晶胞。? 计算设置：合适的3D模型文件一旦确定，必须选择计算类型 和相关参数，例如，对于动力学计算必须确定系综和参数，包括温 度，时间步长和步数。选择运行计算的磁盘并开始CASTAP作业。? 结果分析：计算完成后，相关于CASTAP作业的文档返回用 户，在项目面板适当位置显示。这些文档的一些进一步处理要求获 得可观察量如光学性质。 CASTAP中选择一项任务 1 从模块面板（Module Explorer)选择CASTAP\Calculation。 2 选择设置表。 3 从任务列表中选择所要求的任务。 CASTAP能量任务 CASTAP能量任务允许计算特定体系的总能量以及物理性质。 除了总能量之外，在计算之后还可报告作用于原子上的力；也 能创建电荷密度文件；利用材料观测仪（Material Visualizer)允许目 测电荷密度的立体分布；还能报告计算中使用的Monkhorst-Park的k 点的电子能量，因此在CASTAP分析中可生成态密度图。 对于能够得到可靠结构信息的体系的电子性质的研究，能量任务是 有用的。只要给定应力性质，也可用于计算没有内部自由度的高对 称性体系的状态方程（即压力-体积，能量-体积关系）。注意：具有内部自由度的体系中，利用几何优化（Geometry Optimization)任务可获得 状态方程。CASTAP中能量的默认单位是电子伏特(eV)，各种能量单位的换算 关系见Mohr.P.J(2000).1 eV=0. Ha=23.0605 kcal/mole=96.4853 kJ/moleCASTAP几何优化任务CASTAP几何优化任务允许改善结构的几何，获得稳定结构或 多晶型物。通过一个迭代过程来完成这项任务，迭代过程中调整原 子坐标和晶胞参数使结构的总能量最小化。 CASTAP几何优化是基于减小计算力和应力的数量级，直到小 于规定的收敛误差。也可能给定外部应力张量来对拉应力，压应力 和切应力等作用下的体系行为模型化。在这些情况下反复迭代内部 应力张量直到与所施加的外 部应力相等。 几何优化处理产生的模 型结构与真实结构紧密相似。 利用CASTAP计算的晶格参 数精度列于右图。状态方程计算 在所施加静压力下几何优化可用于确定材料的体模量B和对压力的 导数B‘=dB/dP。过程包括计算理论状态方程（EOS），该方程描述 单胞体积于外部静压力的关系。工艺非常类似于真实实验：使用几 何优化对话框中的应力列表将外部压力固定。通过进行几何优化可 以找到在此压力下的单胞体积。随后的P-V 数据分析与实验研究精 确一致。描述EOS选择分析表达式，其参数适于计算数据点。最流 行的EOS形式是三阶Birch-Murnaghan 方程： P-V 式中V0 为平衡体积。Cohen 等进行了EOS各种解析式的的详细比 较研究。 注意：从相应实验中获得的B和B‘值依赖于计算使用的压力值范围。 利用金刚石压砧获得的实验值通常在0-30GPa范围内，因此推荐理 论研究也在这个范围内。在研究中避免使用负压力值也很重要。此 外，用于生成P-V 数据序列的压力值可能是不均匀的，在低压力范 围要求更精确采样以便获得体模量精确值。几何优化方法 在缺损条件下，CASTAP使用BFGS几何优化方法。该方法通常 提供了寻找最低能量结构的最快途径，这是支持CASTAP单胞优化 的唯一模式。 衰减分子动力学（ Damped molecular dynamics）方法是另一种 可以选择的方法，该方法对具有平滑势能表面的体系如分子晶体或 表面分子与BFGS同样有效。CASTAP动力学任务CASTAP动力学任务允许模拟结构中原子在计算力的影响下将 如何移动。 在进行CASTAP动力学计算以前，可以选择热力学系综和相应 参数，定义模拟时间和模拟温度。 选择热力学系综 对牛顿运动定律积分允许探索体系恒值能量表面（NVE动力学）。 然而，在体系与环境进行热交换条件下发生最本质的现象。使用 NVT系综（或者是确定性的Nosé 系综或者是随机性的Langevin 系综） 可模拟该条件。定义时间步长（timestep ） 在积分算法中重要参数是时间步长。为更好利用计算时间， 应使用大的时间步长。然而，如果时间步长过大，则可导致积分 过程的不稳定和不精确。典型地，这表示为运动常数的系统偏差。注意：量子力学分子动力学计算要求比力场动力学使用更小的时间步长。动力学过程的约束 CASTAP支持Langevin NVT或NVE动力学过程的线性约束。然 而，借助Material Studio界面可以近似使用以下两种更基本的约束： 质心固定 ， 单个原子固定 使用seedname.cell 文档可以利用更复杂的约束。CASTAP性质任务CASTAP性质任务允许在完成能量，几何优化或动力学运行之 后求出电子和结构性质。可以产生的性质如下： ?态密度（DOS）：利用原始模拟中产生的电荷密度和势能，非自 恰计算价带和导带的精细Monkhorst-Pack 网格上的电子本征值。? 能带结构：利用原始模拟中产生的电荷密度和势能，非自恰计算 价带和导带的布里渊区高对称性方向电子本征值。? 光学性质：计算电子能带间转变的矩阵元素。CASTAP分析对话 可用于生成包含可以测得的光学性质的网格和图形文件。 ? 布局数分析：进行Mulliken 分析。计算决定原子电荷的键总数和 角动量（以及自旋极化计算所需的磁矩）。可产生态密度微分计算 所要求的分量。 ? 应力：计算应力张量，并写入seedname.castep 文档。 如果要进行单胞参数固定时进行几何优化运行和要检查点阵偏 离平衡的程度，这些信息是有用的。例如，可进行符合于给定体系 理论基态的固定单胞的点缺陷的超晶胞研究。几何优化后的应力值 显示了与超单胞近似相关联的弹性效应。注意：为计算某种性质，从适当模拟得到的结果文档必须以当前的文件夹形式出现。用第一原理预测AIAs的晶格参数本指南主要是阐明在Materials Studio当中如何 运用量子力学来测定物质的晶体结构。你将从中 学到如何构建晶体结构以及如何设置CASTEP几何 优化运行和分析结果。 本指南的内容如下： 1构建AlAs的晶体结构 2设置和运行CASTEP中的计算 3分析结果 4比较实验数据和结构 注意 如果你的服务器没有足够快的CPU，本指 南限制使用CASTEP进行几何优化计算，因为它会 占用相当长的时间 。1 构建AlAs的晶体结构为了构建晶体结构，我们需要知道你想要构建的晶体的 空间群信息，晶格参数以及它的内部坐标。以AlAs为例， 它的空间群是F-43m或空间群数字是216。它有两种基本 元素Al和As ，其分数坐标分别为（0 0 0）和（0.25 0.25 0.25）。它的晶格参数为5.6622埃。 第一步是构建晶格。在Project explorer的跟目 录上右键单击，选中New | 3D Atomistic Document。 接着在3D Atomistic Document右键单击，把它更名为 AlAs。 从菜单中选择Build | Crystals | Build Crystal，然后显示出Build Crystal对话框，如下：在Enter group中选择 F-43m或在Enter group中单 击，然后键入216，再按下 TAB键.(空间群信息框中 的信息也随着F-43m空间群 的信息而发生变化 )选择Lattice Parameters标签，把a的数 值从10.00改为5.662。单 击Build按钮。一个没有原子的晶格就在3D model document中显示 出来。现在我们就可以添加原子了。 从菜单栏中选择Build | Add Atoms。通过它，我们可以把原子添加 到指定的位置，其对话框如下：在Add Atoms对话框中选择Options标签，确定Coordinate system为Fractional。如上所示。选择Atoms标签，在Element文本框中键入Al，然后按下Add按钮。铝原子就添加到结构中了。 在Element文本框中键入As。在a, b, c文本框中键入0.25。按 Add按钮。关闭对话框。 原子添加完毕，我们再使用对称操作工具来构建晶体结构当中 剩余的原子。这些原子也显示在邻近的单胞中。当然，我们也可以 通过重新建造晶体结构来移去这些原子。 从菜单栏中选择Build | Crystals | Rebuild Crystal...，按下 Rebuild按钮。在显示出的晶体结构中那些原子就被移走了。我们 可以把显示方式变为Ball and Stick。 在模型文档中右键单击，选择Display Styles，按下Ball and stick按钮。关闭对话框。 在3D视窗中的晶体结构是传统的单胞，它显示的是格子的立方 对称。如果存在的话，CASTEP使用的则是格子的全部对称. 既包含有两个原子的原胞和包含有8个原子的单胞是相对应的.不 论单胞如何定义，电荷密度，键长，每一类原子的总体能量都是一 样的，并且由于使用了较少的原子 ,使计算时间得以减少。从菜单栏中选择Build | Symmetry | Primitive Cell。在模型文档中显 示如下： 2设置和运行CASTEP中的计算 从工具栏中选择CASTEP 工具， 再选择Calculation或从菜单栏中选择 Modules | CASTEP |Calculation。CASTEP Calculation对话框如下：AlAs的原胞下面我们将要优化它的几何结构。把Task改为Geometry Optimization ，把Quality改为Fine。 优 化当中的默认设置是优化原子坐标.尽管如此,在本例中我们不仅要 优化原子坐标也要优化晶格.按下Task右侧的More...按钮， 选中Optimize Cell。关闭对话框.当 我们改变Quality时,其他的参数也会有所改变来反映Quality的改变。 选择Properties标签，可从中指定我们想要计算的属性。选中Band structure和Density of states。另外，我们也可以具体指明job control 选项，例如实时更新等。 选择Job Control标签，选中More...按钮。在CASTEP Job Control Options对话框中，把Update 的时间间隔改为30秒。关闭对话框。 按下Run按钮，关闭对话框。 几秒钟之后，在Project Explorer中出现一个新的文件，它包含所有 的运行结果。一个工作日志窗口也会出现，它包含工作的运行状态。我们也可以从Job Explorer中得到工作运行状况的信息。从菜单栏中选择View | Explorers | Job Explorer。Job Explorer中所显示的是与此项目相关联的当正在运行的工作的 状态。它所显示的有用信息有服务器和工作识别数字。如果需要的 话，我们也可以通过Job Explorer来终止当前工作的运行。在工作运行中，会有四个文档打开，它们分程传递关于工作状 态的信息。这些文档包括显示在优化过程中模型更新时的晶体结构， 传递工作设置参数信息和运行信息的状态文档，总体能量图和能量， 力, 应力的收敛以及起重复数作用的位移。工作结束时，这些文件反还给用户。但是由于某个文件可能很大， 也许要耗费教长的时间。 3 分析结果 当我们得到结果是，可能包含数个文档，它们有： (1) AlAs.xsd――最终优化结构。 (2) AlAs Trajectory.xtd――轨线文件，包含每一个优化步骤后的结构。（ (3) AlAs.castep―― 一个输出文本文档，包含优化信息。（ (4) AlAs.param―― 输入模拟信息。对于任一个属性的计算，都会包含有*.param和*.castep文档。运算完毕后，输出结果如下： 其中，第四行的文档为 AlAs.xsd既最终的优化结构。 其 中 第 十 行 为 AlAs.castep 文 档。 在Project Explorer中，单击 AlAs.castep把它击活.在菜单 栏中选择Edit | Find...，在文 本框中键入converged。按下 Find Next按钮。重复按下 Find Next按钮，直到完成搜 索文件的对话框出现。在 Find对话框中，按下Ok和 Cancel按钮。4比较实验数据和结构从我们最初创建的结构单元可以得知，晶格长度应为5.662埃。 所以我们可以比较最小化后的晶格长度和初始时的实验长度。而根 据实验得到的实验长度是基于常规的单元，而非原胞。所以我们需 要我们所创建的单元。 双击AlAs.xsd文件，把它击活。从菜单栏中选择Build | Symmetry | Conventional Cell。 在屏幕上显示出常规的单元。有数种方法可以观测出晶格常数， 但最简单的方法是打开Lattice Parameters对话框。 在模型文档单击右键，选择Lattice Parameters。 其点阵矢量大约为5.629，误差为-0.6%。与实验结果相比较， 这个误差在赝势平面波方法所预期的误差1%―2%之内。在继续下 面的操作之前，先保存项目，再关闭所有窗口。 在菜单栏中选择File | Save Project，然后选择Windows | Close All。 可视化电荷密度 : 从CASTEP Analysis工具中得电荷密度。 从工具栏选择CASTEP ，然后选择Analysis或从菜单栏选择 Modules | CASTEP | Analysis，再选中Electron density选项。此时， 会有一条信息“no results file is available”，所以我们需要明结果文件。在Project Explorer中双击AlAs.castep.。 在Project Explorer中双击AlAs.xsd。从菜单栏中选择Build | Symmetry | Primitive Cell。然后按下Import按钮。结果如下 : 我们可以从Display Style对话框中 改变等能面的设置方式。 在模型文档 中右键单击，选择Display Style，选 中Isosurface标签。Isosurface tab如 下：AlAs的电子密度等能面在这儿，我们可以改变各种设置。 在Iso-value中键入0.1，按下TAB键。注意等能面如何改变。 把Transparency滑块移到右端。当移动透明度滑块是，表面变 的更加透明。 在文档中移动鼠标，旋转模型。当模型旋转时，增加旋转速度 会使等能面以圆点的形式显示出来。从Display Style中，我们也可 以移去等能面。 勾去Visible选项，关闭对话框。 态密度和能带结构 我们可以通过Analysis工具来显示态密度和能带结构的信息。 能带图显示的是在布里渊区中K矢量沿着高对称性方向上的电子能 量依赖度。这些图给我们提供了非常有用的工具，让我们可以对材 料的电子结构进行定性分析――例如，它使我们很容易就可以识别 D态和F态的窄带，与其对立的是近自由电子形成的能带既与S态和P 态相对应的能带。 DOS和PDOS图提供了物质电子结构的优质图象，有时它直接 与实验室的分光镜结果有关。 CASTEP主要的输出文件AlAs.castep，它所包含的能带结构和态密度信息是有的，更加详细的信息包含在AlAs_BandStr.castep文 档中。 打开Analysis对话框，选择Band structure.AlAs_BandStr.castep 文件是自动选上，在次对话框中，我们可以选择在一个图表文件中 同时显示态密度和能带结构。 注意：我们也可以分别分析态密度和能带结构，然后把它们的图形 文档分别显示出来。 在DOS区域，选中Show DOS检验栏。按下View按钮。生成一 个包含态密度和能带结构的图表文件。 我们也可以使用CASTEP来计算许多其它的属性，例如反射率 和非导电性函数等。 ENDCO吸附在Pd（110）面 目的：介绍用CASTEP如何计属表面上的吸附能。模块：CASTEP，Materials Visualizer 背景知识：Pd的表面在许多催化反应中都起着非常重要的作用。 理解催化反应首先是弄清楚分子是如何与这样的表面相结合的。 在本篇文章中，通过提出下列问题，DFT（二维傅立叶变换） 模拟有助于我们的理解：分子趋向于吸附在哪里？可以有多少 分子吸附在表面？吸附能是什么？它们的结构像什么？吸附的 机制是什么？ 我们应当把注意力集中于吸附点,既短桥点，因为众所周知它 是首选的能量活泼点。而且覆盖面也是确定的(1 ML).。在1 ML 覆盖面上CO 分子互相排斥以阻止CO 分子垂直的连接在表面上。 考虑到(1x1)和(2x1)表面的单胞，我们将要计算出这种倾斜对化学 吸收能的能量贡献。 绪论:在本指南中，我们将使用CASTEP来最优化和计算数种系 统的总体能量。一旦我们确定了这些能量，我们就可以计算CO 在Pd(110)面上的化学吸附能。本指南包括：1．准备项目2．最优化Pd 3．构造和优化CO 4．构造Pd(110).面 5 ． RELAXINGPd(110) 面 6 ． 添 加 CO 到 1x1Pd(110) ， 优 化 此 结 构 7．设置和优化2x1Pd(110)面 8．分析能量 9．分析态密度 1．准备项目本指南包含有五种明显不同的计算。为便于管理项目，我们先 在项目中准备五个子文件夹。 在Project Explorer的根图标上右键单击，选择New | Folder。 再重复此操作四次。在New Folder上右键单击，选择Rename，键 入Pd bulk。在其它的文件上重复此操作过程，把它们依次更名为 Pd(110)，CO molecule，, (1x1) CO on Pd(110)，和 (2x1) CO on Pd(110).2．最优化bulk PdMaterials Studio所提供的结构库中包含有Pd的晶体结构。 在Project Explorer中，右键单击Pd bulk文件夹并且选择Import....， 从Structures/metals/pure-metals中导入Pd.msi。 显示出bulk Pd的结构，我们把显示方式改为Ball and Stick。在 Pd 3D Model document中右键单击，选择Display Style，在Atoms 标签中选择Ball and Stick，关闭对话框。 现在使用CASTEP来优化bulk Pd。从工具栏中选择CASTEP ，再选择Calculation或菜单栏 中选择Modules | CASTEP | Calculation。 CASTEP对话框如下：把Task从Energy改为Geometry Optimization，按下More...按钮，在 CASTEP Geometry Optimization对 话框中选中Optimize Cell选项。按 下Run键。出现一个关于转换为原胞 的信息框，按下OK。 工作递交后，开始运行。现在我 们应该进行下一步操作，构造CO分 子。当工作结束后，再返回此处， 显示晶格参数。工作完成后，我们 应保存项目。 选择File | Save Project，然后再从菜单栏选择Window | Close All。 在Project Explorer中打开位于Pd CASTEP GeomOpt文件夹中的 Pd.xsd。显示的即为Pd优化后的结构。 在3D Model document中单击右键，选中Lattice Parameters。其 晶格参数大约为3.91 ?，而其实验植为3.89 ?。3．构造和优化COCASTEP只能处理周期性的体系。为了能够优化CO分子的几何结构，我们必需把它放入晶格点阵中。 在Project Explorer中，右键单击on CO molecule，选择New | 3D Atomistic Document.在3D Atomistic Document.xsd上右键单击，选 中Rename。键入CO，按下RETURN键。 现在显示的是一个空3D模型文档。我们可以使用Build Crystal 工具来创建一个空晶格单元，然后在上面添加CO分子。 从菜单栏中选择Build | Crystals | Build Crystal，再选中Lattice Parameters标签，把每一个单元的长度a, b, 和 c改为8.00，按下 Build按钮。在3D模型文档中显示出一个空单元。 从菜单栏选择Build | Add Atoms。 CO分子中C-O键的键长实验值是1.1283 ?。通过笛卡儿坐标系 来添加原子，我们可以精确的创建此种键长的CO分子。 在Add Atoms对话框中，选择Options标签，确定Coordinate system为Cartesian。然后选中Atoms标签，按下Add按钮。 在Add Atoms对话框中，把Element改为O，x 和 y的坐标值依然 为0，把z的坐标值改为1.1283。按下Add按钮，关闭对话框。 现在我们准备优化CO分子。从工具栏中选择CASTEP工具，然后选择Calculation。先前计算时的设置依然保留着。此次计算不需要优化。 在Setup标签中，按下More...按钮。勾去Optimize Cell选项。关闭 对话框。选择Electronic标签，把k-point set由Medium改为Gamma。 选择Properties标签，选中Density of states。把k-point set改为 Gamma，勾选Calculate PDOS选项。按下Run按钮。 计算开始，我们可以进行下一步操作。4．构造Pd(110).面下面我们将要用到从Pd bulk中获得的Pd优化结构。 从菜单栏中选择File | Save Project，然后在选中Window | Close All。在Pd bulk/Pd CASTEPGeomOpt文档中打开Pd.xsd。 创建表面分为两个步骤。第一步是劈开表面，第二步是创建一 个包含表面的真空板。 从菜单栏中选择Build | Surfaces | Cleave Surface，把the Cleave plane (h k l)从（-1 0 0）改为（1 1 0），然后按下TAB键。把 Fractional Depth增加到1.5，按下Cleave按钮，关闭对话框。 此时，显示出一个包含有二维周期性表面的全新的三维模型文档。尽管如此，CASTEP要求有一个三维周期性的输入体系。我们可以 Vacuum Slab工具来获得。 在菜单栏中选择Build | Crystals | Vacuum Slab，把Vacuum thickness从10.00改为8.00。按下Build键。 则结构由二维变成三维，把真空添加到了原子上。在继续下面的 操作前，我们要重新定位一下格子。我们应该改变格子的显示方式 并且旋转该结构，使屏幕上的Z轴成竖直状。 在3D model document中单击右键，选择Lattice Parameters选项。 选择Advanced标签，按下Reorient to standard按钮，关闭对话框。 在3D model document中单击右键，选择在 Display Style。然后选中Lattice标签，在Display中， 把Style从Default改为Origina。 按下Up 指针键两次，三维模型文档如右所示： 把Z坐标最大值所对应的Pd原子称为最高层Pd原子。 在本指南的稍后部分，我们要求知道原子层间 的距离do，我们可以通过计算原子坐标来得到。 从菜单栏中选择View | Explorers | Properties Explorer，选择FractionalXYZ中X=0.5，Y=0.5 的 Pd原子。注意从XYZ属性中所获得的Z的坐标值。Z的坐标值应为1.386 ?，此既为原子层间的距离。 注意：一个fcc(110)体系，do 可通过下列公式得到： . 在释放表面之前，如果仅仅是只需要释放表面，我们必需要束 缚住Pd原子。 按住SHIFT键选中所有的Pd原子，不包括最高层的Pd原子。从 菜单栏中选中Modify | Constraints，勾选上Fix fractional position。 关闭对话框。则刚才所选中的原子已经被束缚，我们可以通过改变 显示的颜色来看到它们。 在3D模型文档中单击以取消所选中的原子。右 键单击选择Display Style，在Atoms标签的Coloring 部分，把Color by选项改为Constraint。3D模型文 档显示如右： 把Color by选项再改为Element，关闭对话框。 从菜单栏中选择File | Save As...，把它导引到Pd(110) 文件夹中，按下Save按钮。对(1x1) CO on Pd(110)文 件夹也重复此操作，但是这一次把文档的名字改为 (1x1) CO on Pd(110)。再选择File | Save Project，然 后再选择Window | Close All。5．释放Pd(110)面现在我们最优化Pd (110)表面。 在Project Explorer的Pd (110)文件夹中打开Pd(110).xsd。从工具 栏中选中CASTEP 工具，然后选择Calculation。按下More... 按钮，确定Optimize Cell没有被选中。关闭对话框。 为了维持我们想要完成的计算的连贯性，我们应该更改Electronic 标签中的一些设置。 选 择 Electronic tab 标 签 ， 然 后 按 下 More... 按 钮 。 从 CASTEP Electronic Options 对 话 框 中 选 择 Basis 标 签 ， 勾 选 上 Use custom energy cut-off并且把域植从260.0改为300.0。选择k-points标签，勾 选上Custom grid参数。在Mesh parameters域中，把a改到3，b改到 4，c改到1。关闭对话框。 我们还应该计算此体系的态密度。 选择CASTEP Calculation对话框中的Properties标签，选中 Density of states。勾选上Calculate PDOS，把k-point set改为 Medium。 按下Run按钮，关闭对话框。 计算的运行会耗费一定的时间，我们可以最后做分析。我们现在可以构建下一组表面。 从菜单栏中选择File | Save Project，选择Window | Close All。6．添加CO到1x1Pd(110)，优化此结构我们要使用在(1x1) Co on Pd(110)文件中的结构来进行下面的工作。 在Project Explorer中，打开(1x1) Co on Pd(110)文件中的(1x1) CO on Pd(110).xsd。 现在在short bridge position上添加CO分子。我们要利用的依据是： CO 在 Pd(110)上的键长已经通过实验所获得。右图中阴影线原子在格子中不显示： Original display mode。 第一步是添加碳原子。Pd-C键的键长 （用dPd-C表示）应为1.93 ?。当我们使 用Add Atom 工具时，我们即可以使用 笛卡儿坐标也可以使用分数坐标，但在 本例当中，我们应该使用分数坐标xC, yC, 和zC。xC, yC非常简单，xC =0，yC =0.5。 尽管如此，zC比较困难。我们可以通过 zPd-C 和zPd-Pd二者之间的距离来构造它。CO在Pd(110)的yz平面上 的几何结构zPd-Pd可以由晶格参数a0除以√2得到（它应为2.77 ?）。zPd-C可从公式得到（它应为1.35 ?）。把zPd-C 和 zPd-Pd相加可获得zC（它应为4.12 ?）。现在我们把距离改 为分数长度，可以通过晶格参数（Lattice parameters）工具得到。 在3D模型文档中单击右键，选择Lattice parameters。注意c的值。 为了计算z的分数坐标，我们仅需要用晶格参数c除以zC（结果为 0.382 ?）。 从菜单栏中选择Build | Add Atoms，然后选中Options标签。确保 Coordinate system为Fractional。选择Atoms标签，把a改为0.0，b 为0.5，c为0.382。按下Add按钮。 如果我们想确认我们已经正确的设置了模型，可以使用 Measure/Change工具。 单击工具栏中Measure/Change工具 的选项箭头，然后选择 Distance。在Pd-C键上单击。 下一步是添加氧原子。 在Add Atoms对话框中，把Element改为O。 在实验中，C-O键的长度为1.15 ?。在分数坐标中它为0.107，把 这个值添加到碳的z分数坐标上（0.382），氧的z坐标值为0.489。把c的域值改为0.489，按下Add按钮。关闭对话框。 Pd最原始的对称性是P1，但随着CO的添加它以改变。我们可以 通过Find Symmetry工具来找到其对称性和强加对称性（Impose Symmetry）。 在工具栏中选择Find Symmetry工具 ，按下Find Symmetry按钮，随后按下Impose Symmetry按钮。对称性为PMM2。 在3D模型文档中单击右键，选择Display Style。 选中Lattice标签，把Style改为Default。 结构如右所示： 在优化几何结构之前，我们先把它保存到(2x1) CO on Pd(110)文件夹中。 从菜单栏中选择File | Save As...，引导到(2x1) CO on Pd(110)文件。把文档保存为(2x1) CO on Pd(110).xsd。现在可以优化结构 。 从菜单栏中选择File | Save Project ，然后选择 Window | Close All。在Project Explorer中，打开 (1x1)CO on Pd(110)文件夹中的(1x1)CO on Pd(110) .xsd。从工具栏中选择CASTEP 工具， 然后选择Calculation。从先前的计算中得到的参数应当保留。 按下Run按钮。 在运行过程中，我们进行最后结构的构建。7．设置和优化2x1Pd(110)面第一步是打开(2x1) CO on Pd(110)文件夹中的3D模型文档。 在Project Explorer中，打开(2x1) CO on Pd(110)文件夹中的(2x1) CO on Pd(110).xsd。 这就是当前的1x1单元，我们需要使用Supercell工具把其变为 2x1单元。 从菜单栏中选择Build | Symmetry | SuperCell， 把b增加到2，按下Create Supercell按钮。关闭对 话框。其结构看起来如右： 现在我们使CO分子倾斜。 为了简化此操作, 定义位于 y=0.5处的分子为 A分子，位于y=0.0 处的分子为B分子。 选择B分子的碳原子。在Properties Explorer 中，打开XYZ属性，从X域中减去0.6。 对于B分子的氧原子重复此操作，但从X域 (2x1) Cell of CO on 中减去1.2。 Pd(110)选择A分子的碳原子。在Properties Explorer中，打开XYZ属性， 在X域中增加0.6。对于A分子的氧原子重复此操作，但在X域中增 加1.2。分子的X轴向下的视图如右： 尽管如此，我们应注意到Pd-C和CO键长的最初值已经改变。 选中A分子的碳原子，使用Properties Explorer，把FractionalXYZ属性中Z的 域值改为0.369。对B分子重复上述操作。 此操作在于更正Pd-C的键长。我们可以 使用Measure/Change工具来更正C-O键长。 在工具栏中单击Measure/Change工具 的选项箭头，选中 Distance。单击A分子的C-O键，在工具栏中选择3D Viewer Selection Mode工具 ，选择监视窗口（既3D Atomistic Document）。在Properties Explorer中，改变Filter to Distance。把 Distance属性改为1.15 ?。对B分子重复此操作。 现在重新计算此体系的对称性。 在工具栏中选择Find Symmetry工具，按下Find Symmetry按钮， 随后再按下Impose Symmetry按钮。现在它的对称性是PMA2。下面我们来优化它的几何结构。 从工具栏中选择CASTEP工具，然后选择Calculation。 对于本次计算，我们需要改变k点的格子参数，这样我们可以比较 本次计算和上次计算的能量值。 选中CASTEP Calculation对话框中的Electronic标签，按下More... 按钮。选择k-points标签，把Custom grid parameters改为：a = 2, b = 3, c = 1。关闭对话框，按下Run按钮。 计算开始。计算结束后，在下面的内容中我们需要详细的摘录整个 体系的能量。我们可以进行下一步，摘录先前计算的能量。8．分析能量在这一部分，我们将要计算化学吸收能DEchem，定义如下：允许CO分子依着彼此倾斜，然后减低分子的自我排斥力，会导致能 量的增加。排斥能可从下面的公式得到： 为计算这些属性，我们需要从CASTEP的文本输出文档中摘录每一 次模拟的整个能量。在Project Explorer中，打开CO molecule/CO CASTEP GeomOpt文件 夹中的CO.cst。按下CTRL+F键，搜索Final electronic.。向下滑动数 行，在下面的表格中记录下出现在“TOTAL ENERGY IS”此行之后 的数值。重复此操作，找到其它体系的整个能量，完成下面的表格。如果获取了所有的能量值，使用上面的等式很简单便可计算出 DEchem and DErep。他们的数值大约分别为1.9eV和72MeV。9．分析态密度下面我们要检查态密度（DOS）的改变。这会使我们对CO在 Pd(110).上的连接机制有更深入的了解。为了做到这一点，我们需 要显示孤立的CO分子和(2x1) CO 在Pd(110)上的态密度。 在Project Explorer中，打开CO molecule/CO CASTEP GeomOpt文 件夹中的CO.xsd。 从工具栏中选择CASTEP工具，然后选择Analysis，选中Density of states，选上Partial，不选f 和sum，但其他的选项都保持原先的状 态。按下View按钮 。显示出CO分子PDOS的图表文档。如右：对(2x1) CO on Pd(110).xsd重复上面 的操作。PDOS of CO moleculePDOS of (2x1) CO on Pd(110) 很明显，孤立的CO分子的电子态大约在-20eV, -5eV 和 -2.5eV， 比CO约束在表面时的能量要低一些。背景 探索任一反应的势能面需要反应过程中每一步的结构和能量（或动 力学和热力学）的快照。尤其重要的是决定反应速度的步骤，它常 常涉及到决定着另人难以捉摸的过度态结构。有许多技术被用来寻 找过度态结构，其中非常出名和有效的是（Linear Synchronous Transit ）LST 和（Quadratic Synchronous Transit）QST方法。 本指含盖内容如下： 1 设置计算的结构 2 优化几何结构 3 定义原子配对 4 使用LST/QST/CG方法计算过度态使用CASTAPLST/QST工具进行度态搜索1 设置计算的结构构建Pd (1 1 1)表面，先导入Pd晶体结构。 从菜单栏中选择File | Import。在structures/metals/pure-metals中选 择Pd.msi。 现在更改此结构的显示方式。 在Pd.xsd中单击右键，选中Display Style。从Atom标签的显示方 式列表中选择Ball and Stick。关闭对话框。Materials Visualizer中的Cleave Surface工具允许我们劈开任一个大块 晶体的表面。 从菜单栏中选择Build | Surfaces | Cleave Surface。 Cleave Surface的 对话框如右： 把Cleave plane中的米勒指数从（-1 0 0）改为（1 1 1）。把Fractional Depth设 置为2.0。按下Cleave按钮。 选择Surface Mesh标签，设置表面向 量U为0.5 -1 0.5，然后按下TAB键。再 设置表面向量V 为0.5 0.5 -1，然后按下 TAB键。关闭对话框。 于是打开了一个包含2D周期性表面的新 的3D模型文档。尽管如此，CASTEP需 要一个作为输入的3D周期性体系。我们 可以使用Vacuum Slab工具获得它。 从菜单栏中选择Build | Crystals | Build Vacuum Slab，把Vacuum thickness从10.00改为7.00。按下Build按钮。 此结构从二唯周期性变成三唯周期性结构，在原子上添加了一层真 空。我们可以移去单胞底部的对称性图形 ，对称性图形同时也出现在晶胞的上部。 从菜单栏中选择Build | Bonds，在Bonding Scheme标签中勾选上 Monitor bonding。关闭对话框。在Pd (1 1 1).xsd中右键单击，选中 Display Style，选择Lattice标签，把Style 设置为In Cell。最后，再把 Style 设置为Default。关闭对话框。 现在我们可以使用已经建造好的Pd (1 1 1)面去构造与反映物所对应 的结构。Pd (1 1 1)面如下所示： 从菜单栏中选中File | New，选择 3D Atomistic Document。当出现提示 时，保存对Pd (1 1 1).xsd所坐的修改。 一个新的空文档出现。选中Pd (1 1 1).xsd把它激活。从菜单栏中选择 Edit | Select All，接着再选择Edit | Copy。 在Project Explorer中选择3D Atomistic Document.xsd把它激活。从 菜单栏中选择Edit | Paste。Pd (1 1 1)晶 体结构出现在新文档中。在文档中某 处单击一下取消所选中的图形。在Project Explorer的3D Atomistic Document.xsd上右键单击，选择 Rename，键入reactants。 添加氢原子： 从菜单栏中选择Build | Add Atoms。 使用Add Atoms工具可以把 原子添加在晶胞指定的位置。 Add Atoms对话框如右： 在 Element 文 本 框 中 ， 键 入 H 。 选 择 Options标签，把Coordinate System设置为 Fractional。返回到Atoms标签，设置a 为 0.56, b 为0.47 和c 为0.70。按下Add按钮。 一个氢原子出现在晶胞中 。 使用相同的步骤，把第二个氢原子添 加到a = 0.47, b = 0.56和c = 0.70位置。关闭 对话框。 提示：当我们添加第二个氢原子时， Materials Studio会产生一个警告信息。这个 警告信息之所以会出现是因为我们所添加的第二个氢原子在第一个 氢原子所定义的公差范围内。在这种情况下，我们建造一个H-H键 长小于1.0?的氢分子既可。选择Yes继续添加氢原子 。一个H2分子是由键长为0.743?的H-H键所形成的。氢分子位于和PdPd键平行的晶胞中心，距离表面大约4.00?。 为了简单起见，我们假设在反应期间表面是固定的。为了做到这一 点，我们必须约束表面原子保留在当前位置。 选择reactants.xsd中的一个Pd原子，然后按下ALT键，再双击选 中所有的Pd原子。从菜单栏中选择Modify | Constraints，确定Pd原子 的坐标系为分数坐标或笛卡儿坐标 。 现在我们来建造产物的结构。这次，我 们要以reactants.xsd的结构为起点。 从菜单栏中选择File | New，再选中 3D Atomistic Document。一个空文档出 现。选择reactants.xsd将其激活。从菜 单栏中选择Edit | Select All，随后再选 择Edit | Copy。 在Project Explorer中选择3D Atomistic Document.xsd将激活。从菜单 栏中选择Edit | Paste。反应物结构出现 在文档中。在文档中某处单击以取消所 反应物结构 选中的图形。在Project Explorer中的3D Atomistic Document.xsd中右键单击， 选择Rename，键入products。 建造产物 ： 在这一部分我们要使用Properties Explorer来改变结构中氢原子的位 置。 从菜单栏中选择View | Explorers | Properties Explorer。 在products.xsd中的一个氢原子上单击。在Properties Explorer中， 显示处FractionalXYZ坐标 。查找位于0.47 0.56 0.70的氢原子，在 FractionalXYZ的文本框中双击，把分数坐标值改为0.67 0.414103。通过相同的步骤把位于0.56 0.47 0.70的氢原子移动到 0...414103。 从菜单栏中选中Build | Bonds，勾选 上Bonding Scheme标签中的Monitor bonding。关闭对话框。 在此新结构中，晶胞中心的两个Pd原 子每一个都有一个氢原子与其连接， 距离大约为1.583?。其图形如右：注意：反应物和产物具有相同的晶格参数。这是必须的，既然在 CASTEP中应用的LST/QST版本不考虑晶格参数的改变。尽管如此， 只要你感兴趣的是那些晶胞改变并不重要的过程，例如在表面的反 应，原子扩散和在体积中的空缺等，那么它的局限性就不是很明显。2 优化几何结构正如我们已经讨论的那样，反应物和产物的cell parameters必须相同。 基于此种原因，任何优化仅涉及晶胞种原子的位置。而且我们可以 认为反应物的结构和原先一样，所以我们不需要优化它们的结构。 确定products.xsd文档处于激活状态。从工具栏中选择CASTEP 工具 ， 然 后 选 择 Calculation 或 者 从 菜 单 栏 中 选 择 Modules | CASTEP | Calculation。 CASTEP Calculation对话框如右所示： 下面我们开始优化它的几何结构。 把Task 改为Geometry Optimization。把 Quality改为Medium。 选择Electronic标签，把k-point set设置 为Gamma。 按下More...按钮，然后选择SCF标签。 把Charge设置为0.4。按下和此选项先相关的More...按钮，把DIIS history list 改为5。 返回SCF标签，确定没有选上Fix occupancy选项。关闭对话框。 我们也可以指定工作控制选项，例如实时更新。 选择Job Control标签。按下More...按钮，在CASTEP Job Control Options对话框中，把Update interval 改为30.0秒。关闭对话框。 如果你在另外的计算机上运行工作，你同样可以在Job Control标签 中如此选择。 按下Run按钮。关闭CASTEP Calculation对话框。 很快，在Project Explorer中出现了一个新文档。它包括计算的所有 结果。一个Job Log窗口显示出来，它包括工作的状态。你也可以从 Job Explorer中获得此信息。 从菜单栏中选择View | Explorers | Job Explorer。 Job Explorer显示的和当前项目相关的任何激活的工作的状态。它所 显示的有用信息包括服务器和工作识别数字。我们也可以使用Job Explorer来停止工作。 当工作进行时，打开了四个关于工作状态的文档，它们分程传递信 息。这些文档包括显示在优化过程中模型更新时的晶体结构，传递 工作设置参数信息和运行信息的状态文档，总体能量图和能量， Forces, Stress 的收敛以及起重复数作用的位移。3 定义原子配对对CASTEP来说，为了完成过度态搜索，反应物文档和产物文档中 的所有原子都需要配对。此任务可以使用Reaction Preview工具来完 成，此工具可从工具栏中得到。 第一步，我们应并排显示反应物和产物。 从菜单栏中选择Window | Tile Vertically。 现在，我们开始使反应物和产物中的原子配对。 从菜单栏中选择Tools | Reaction Preview。 其对话框如右： 分别选择reactants.xsd和products.xsd为 反应物和产物。单击Match...按钮。 显示出Find Equivalent Atoms对话框。 我们应从最开始看此对话框，没有原子 匹配和8个原子不能匹配。现在有6个匹 配原子和两个不能匹配原子。这两个不 能匹配原子为氢原子。 在反应物列表中双击2xH。 反应物栏中与之对应的文件夹被打开。 反应物栏应包括7:H和8:H。反应物和产物的3D模型文档都打开后，在反应物窗口中单击 7:H。在产物窗口中单击7:H。两个窗格中的氢原子应该都被选中， 这些氢原子在两个3D模型文档中应该使一样的。 如果它们使一样的，单击Set Match。 我们可以回想以下反应物和产物的匹配。 单击反应物栏或产物栏列表中的一个原子，回顾一下其恰当的 配对。如果对它们的匹配满意的话，关闭对话框。 为了利用CASTEP LST/QST完成过度态搜索，我们需要创建反映物 和产物之间的路径就如同CASTEP计算时的输入一样。 在Reaction Preview对话框中，确定Number Of Frames是10。按 下Preview按钮。关闭Reaction Preview对话框。 数秒钟之后，显示出一个名为reactants-products.xtd,的新的3D轨线文 档。我们可以对这个文件运行CASTEP计算。通过使用 Animation工 具栏我们可以使轨线文档充满生气。 如果看不到的话，使用View菜单显示出Animation toolbar。使用 Play按钮使轨线充满生气。在Bounce模式下，效果最好。 我们可以转换为键监视模式，这样的话每一步之后，键都会重新计 算。 从菜单栏中选择Build | Bonds，勾选上Bond Monitoring。关闭对话 框。再一次使轨线充满生气。当看完之后，按下Stop按钮。4 使用LST/QST/CG方法计算过度态现在我们设置CASTEP计算以确定过度态。 从Modules工具栏中大开CASTEP Calculation对话框。在Setup标 签中，把Task从Geometry Optimization改为TS Search。 单击More...按钮，确定Search Protocol设置为Complete LST/QST 以及Quality为Medium。关闭对话框。 电子的Hamiltonian函数设置和几何优化计算时的一样。最后，我们 应更新工作描述。 选择Job Control标签，确定没选上Automatic检验栏。键入TS作 为Job description。 按下Run按钮，关闭CASTEP Calculation对话框。 在计算过程中，在Visualizer里显示出几个不同的文档和一个 LST/QST图形。这些是报告计算的状态。特别是LST/QST图通过显 示一张能量经由LST, QST and CG的路径坐标图来监视TS搜索过程。 过渡态计算完成后，过渡态写在TS.xsd文档里。 如果文档没有自动显示出来，在Project Explorer中双击TS.xsd 。 计算的文本结果包含在TS.castep文件中。如果文档没有自动的显示出来，双 击Project Explorer中的TS.castep。按下 CTRL+F搜索Energy barrier。 反应物的能量大约是-0.96133 eV。从 反应物到过渡态的势垒大约为2.97645 eV，从过渡态到产物的势垒大约为3.93778 eV。过渡态结构计算BN的弹性常数背景 当前可应用于大周期性体系的密度功能理论（DFT）方法的发 展已经成为解决材料设计和处理难题的重点。这个理论可以使我们 阐释实验数据和预测材料潜在的属性，从一种未知晶体的结构属性 到结合能和表面分子的反应都可以预测。这些工具可以用来指导和 引导我们设计新材料，允许研究人员推定潜在的化学和物理过程。 绪论 在本指南中，我们将学习如何使用CASTEP来计算弹性常数和 其他的力学性能。首先我们要优化BN立方晶体的结构，然后计算它 的弹性常数。 本指南主要包括以下内容：1 优化BN立方晶体的结构2 计算BN的弹性常数 3 弹性常数文件的描述1 优化BN立方晶体的结构在计算弹性常数之前并不是必须要进行几何优化，所以对于实 验观测到的结构我们可以产生Cij数据。尽管如此，如果我们完成几 何优化包括晶胞优化，可以获得更多相容的结果，然后再计算与理 论基态对应的结构的弹性常数。 弹性常数的精确度，尤其是切变常数的精确度，主要取决于 SCF计算的品质，特别是布里渊区取样和波函数收敛程度的品质。 所以我们应该使用SCF公差为Fine的设置和 k点取样以及来自FFT格 子的Fine设置。 首先应导入BN结构。 在菜单栏中选择File | Import， 从 structures/semiconductors 中 选中BN.xsd。BN的晶体结构如右： 现在设置几何优化。 从工具栏中选择CASTEP 工 具，然后选择Calculation或从菜单栏 中选择Modules | CASTEP | Calculation。 CASTEP Calculation对话框如下：在Setup标签中，把Task设置为 Geometry Optimization，把Quality 设置 为Fine，并且把Functional设置为GGA and PW91。 选择Electronic标签，按下More... 按 钮 以 得 到 CASTEP Electronic Options对话框。把Derived grid的设 置从Standard改为Fine。关闭CASTEP Electronic Options对话框。 选择Job Control标签，选中你想要 运行CASTEP工作的Gateway。 按下CASTEP Calculation对话框中的Run按钮。 优化之后，此结构的晶胞参数应为a=b=c=3.60556?。现在我们可以 继续计算优化结构的弹性常数。2 计算BN的弹性常数选择CASTEP Calculation对话框中的Setup标签。从Task的下拉清 单中选择Elastic Constants。按下More...按钮。 CASTEP Elastic Constants对话框如下：把Number of steps for each strain从4 增加到6。关闭CASTEP Elastic Constants对话框。按下CASTEP Calculation对话框中的Run按钮。 CASTEP的弹性常数计算任务的结果以 一批.castep输出文件的形式给出。这 些文件中的每一个文件都代表确定的 晶胞在假设的应变模式和应变振幅下 的几何优化运行结果。对这些文件的 习 惯 命 名 方 式 为 ： seedname_cij__m__n。对于给定的模 式来说，m代表当前的应变模式，n代 表当前的应变振幅。 CASTEP可以使用这些结果来分析 每一个运行计算出来的压力张量和产生一个有关弹性属性信息的文 件。 从工具栏中选择CASTEP 工具，然后选择Analysis或者 从菜单栏中选择Modules | CASTEP | Analysis。 从属性清单中选择Elastic constants，从BN的弹性常数计算工作中得到的结果文件应自动显示在Results file选框中。按下Calculate按 钮。 在结果文件夹中创建了一个新的文档BN Elastic Constants.txt。 此文档中的信息包括输入的应变摘要和计算出的应力；每一种 应变模式的线性适配结果包括适配的质量；在计算出的应力和给定 对称性的弹性常数之间的对应;弹性常数Cij和弹性柔量Sij的表格以及 最后推知的属性如体积模量和于之相反的可压缩性杨氏模量以及三 唯方向上的Poisson比率和建造一个物质为各向同性媒介的模型所需 要的Lame 常数。3 弹性常数文件的描述对于这种点阵类型需要两种应变模式。对于每一种应变模式， 这里都有一个从各自的.castep文件中获取的计算出的压力的摘要。=============================================== Elastic constants from Materials Studio: CASTEP =============================================== Summary of the calculated stresses ********************************** Strain pattern: 1 ====================== Current amplitude: 1 Transformed stress tensor (GPa) : -12.........365274 Current amplitude: 2 Transformed stress tensor (GPa) : -13.........545781提供了应力，应变的组成和弹性常数张量之间联系的所有信息。在 这一阶段，每一个弹性常数均有一个简洁的指数代表而不是由一对 ij指数代表。稍后会在文件夹中给出压缩符和常规的指数标定之间的对应。 和弹性系数相对应的应力 (压缩符)： 1 7 7 4 0 0 as induced by the strain components: 1 1 1 4 0 0 在下面的表格中给出了每一种应力组成的应力-应变线性适配关系:Stress Cij value of value of index index stress strain 1 1 -12.. 1 -13.. 1 -14.. 1 -15.. 1 -16.. 1 -17..003000 C (gradient) : 768.968762 Error on C : 5.954823 Correlation coeff: 0.999880 Stress intercept : -14.807379此梯度提供了弹性常数的数值（或弹性常数的线性组合），适配的 质量，由相关系数表示，提供了另人满意的弹性常数的不确定度。在进一步的分析中没有使用压力的切点值, 它很简单的指示出收敛 的基态离最初的结构有多远。 所有应变模式的结果总结如下：============================ Summary of elastic constants ============================ id i j Cij (GPa) 1 1 1 768.96876 +/- 5.955 4 4 4 443.89917 +/- 1.004 7 1 2 142.23838 +/- 1.569The errors are only provided when more than two values for the strain amplitude were used, since there is no statistical uncertainty associated with fitting a straight line to only two points.弹性常数以常规的6x6张量的形式显示出,随后弹性柔量 （compliances）以相似的6x6形式显示出：===================================== Elastic Stiffness Constants Cij (GPa) ===================================== 768.838 142.00 0.00 142.876 142.00 0.00 142.838 768.00 0.00 0.00 0.917 0.00 0.00 0.00 443.00 0.00 0.00 0.917 ======================================== Elastic Compliance Constants Sij (1/GPa) ======================================== 0....................................0022528文件的最后部分包含推出的属性:Bulk modulus = 351.14851 +/- 2.244 (GPa) Compressibility = 0.00285 (1/GPa) Axis Young Modulus Poisson Ratios (GPa) X 724.56227 Exy= 0.1561 Exz= 0.1561 Y 724.56227 Eyx= 0.1561 Eyz= 0.1561 Z 724.56227 Ezx= 0.1561 Ezy= 0.1561Lame constants for isotropic material (GPa) Lambda = -118.8296, Mu = 443.8992END预测锗的热力学属性背景 线性响应或密度功能混乱理论是点阵动力学从头开始计算中最 受欢迎的方法之一，尽管如此，这种方法的应用已经扩充到对振动 属性的研究。线性响应提供了一种分析方法用于计算给定混乱的二 级派生的整体能量。可以计算出许多属性，主要依赖于混乱的种类。 在离子位置的混乱可以引起动力矩阵和声子；在磁场中引起NMR效 应；在单位晶格矢量中产生弹性常数；在电场中引起非传导性效应 等。 在本指南中，我们将要学习为了计算声子散射和能态密度以及 预测热力学属性如焓和自由能，如何使用CASTEP来完成线性响应 计算。 本指南主要包含以下内容：1 优化锗单胞的结构2 计算声子散射和能态密度 3显示声子散射和能态密度 4显示热力学属性1 优化锗单胞的结构首先我们要导入锗的结构，它包含在Materials Studio所提供的结 构库中。 在菜单栏中选择File | Import。遵循下列路径 structures/metals/pure metals选中Ge.xsd。 把它转换为原胞结构后，对它的计算会更快。 从 菜 单 栏 中 选 择 Build | Symmetry | Primitive Cell。锗的原胞结构如右所示： 现在开始使用CASTEP来优化锗的几何结 构。 从工具栏中选择CASTEP工具 然 后 选 中 Calculation 或 从 菜 单 栏 中 选 择 Modules | CASTEP | Calculation 。 CASTEP Calculation的对话框如下：几何优化的默认值不包括对单胞的优 化。 在Setup标签上，把Task从Energy 改为Geometry Optimization，把 Functional改为LDA。 在CASTEP Geometry Optimization 对 话框中，按下More...按钮，勾选上 Optimize Cell。 选 中 Electronic 标 签 ， 把 Energy cutoff 设 置 为 Ultra-fine ， 把 SCF tolerance设置为 Ultra-fine，把kpoint set 设 置 为 Coarse 以 及 把 Pseudo-potentials设置为 Norm-conserving。 选中Job Control标签。选择你想要在其上运行工作的Gateway location。把Runtime optimization设置为Speed。 按下Run按钮开始运行。 工作递交后开始运行。它大概需要2分钟时间，这主要取决于你的 电脑的速度。结果放在一个被称为Ge CASTEP GeomOpt的新文件夹 中。2 计算声子散射和能态密度为了计算声子散射和声子的能态密度，在从CASTEP Calculation对 话框的Properties标签选定适当的属性后，我们必须完成一个单点 能量计算。 确定Ge CASTEP GeomOpt文件夹中的Ge.xsd文件示激活文档。 选中CASTEP Calculation对话框中的Setup标签，把Task设置为 Energy。 选 择 Electronic 标 签 ， 按 下 More… 按 钮 ， 显 示 出 CASTEP Electronic Options对话框。选择SCF标签，勾选上Fix occupancy。 关闭CASTEP Electronic Options对话框。在CASTEP Calculation 对话框中选定Properties标签。选择Phonon dispersion并且把qvector set设置为medium。选择Phonon density of states并且把 q-vector set设置为medium。 选 中 Job Control 标 签 。 选 择 你 想 要 在 其 上 运 行 工 作 的 Gateway location。 按下Run按钮。 在Ge CASTEP GeomOpt文件夹中创建了一个名为Ge CASTEP Energy的新文件。当能量计算完成后，两个新文件 Ge_PhononDisp.castep和Ge_PhononDOS.castep放在此文件夹中。3显示声子散射和能态密度声子散射曲线显示出声子能量沿着布里渊区高对称性方向如何依赖 于q向量。此信息可以从单晶的中子散射实验中获得。只有为数不 多的物质可以获得这样的信息，所以用来确定建模方法是否正确的 理论偏差曲线对于论证在卷首的从头开始计算方法的预测性能力是 非常有用的。在一定情形下，它可能测量态密度而不是声子散射。 而且与声子的态密度有直接关系的电子―声子交感作用可以通过隧 道实验直接测量。所以能够从第一原理计算出声子的态密度是非常 重要的。Materials Studio可以从任何.phonon CASTEP输出文件中 产生声子散射图和态密度图。这些文件隐藏在Project Explorer里， 但 是 一 个 .phonon 文 件 会 和 每 一 个 拥 有 PhonDisp 或 PhonDOS 后 缀 的.castep文件一起产生。 现在，我们使用先前的计算结果来创建声子散射图。 从 Materials Studio 的 菜 单 栏 中 选 择 Modules | CASTEP | Analysis。从属性列表中选择 Phonon dispersion。确定Results file选择框中显示的是Ge_PhononDisp.castep。 从Units下拉列表中选择cm-1。从Graph style下拉列表中选择Line。 按下View按钮。在 结 果 文 档 中 创 建 了 一 个 新 的 图 形 文 档 Ge Phonon Dispersion.xcd。它应和下图相似：声子散射的实验图如下所示：预测的频率可从Ge_PhononDisp.castep文件中得到。 在Project Explorer中双击Ge_PhononDisp.castep。按下CRTL+F 键，搜索Vibrational Frequencies。 结果文件中的部分内容如下所示：======================================================= + Vibrational Frequencies + + -------------------------------+ + + + Branch number Frequency (cm-1) + + ==================================================== + + + + ------------------------------------------------------------------------------------------- + + q-pt= 1 ( 0....076923 + + ------------------------------------------------------------------------------------------- + + 1 122.413706 + + 2 122.413706 + + 3 211.138305 + + 4 211.138305 + + 5 288.235680 + + 6 288.235680 + + ------------------------------------------------------------------------------------------- + + q-pt= 2 ( 0....038462 + + -------------------------------------------------------------------------------------------- + ....每一个q点和每一个分支（纵向光波或声波(LO/LA),横向光波或声 波(TO/TA)的频率以cm-1 给出，同时也给出了q点在倒易空间中的位 置。高对称性点G, L和X在倒易空间中的位置各自为(0 0 0), (0.5 0.5 0.5) 和 (0.5 0 0.5)。这些点和q点12, 6 以及19相对应。 预测的频率（cm-1）和实验的频率（cm-1）如下：总体来说，计算的精度是可以接受的。在Gamma点错误的声学波频 率，3 cm-1 而不是0cm-1，使我们感觉这只是一般的精确度。通过对 更好的SCF k点格子的计算，我们可以获得更加另人满意的实验结 果。 现在创建声子态密度图： 从Materials Studio 菜单栏中选择Modules | CASTEP | Analysis 。 从属性列表中选择Phonon density of states。确定Results file 选择框 中显示的是Ge_PhononDOS.castep。把DOS display设置为Full。按下 View 按钮。 创建了一个新的图形文档Ge DOS.xcd。它应当与下图相似：4显示热力学属性在CASTEP中的声子计算可以用来评价近似准谐波晶体的焓，熵， 自由能，格子的热容对于温度的依赖性。可以用这些结果和实验数 据（如，热容的测量）相比较以预测不同的结构经过修正后的相稳 定性和相转变。 所有与能量相关的属性均画在同一种曲线图中，并且0点能量 的计算值也包括在内。热容被独自画在图表的右侧。 现在使用声子计算的结果创建热力学属性图表。 从Materials Studio 菜单栏中选择Modules | CASTEP | Analysis 。 从属性列表中选择Thermodynamic properties。确定Results file 选择 框中显示的是Ge_PhononDOS.castep。 勾选上Debye temperature图，按下View按钮。 在结果文件夹中创建了两个新的图形文档Ge Thermodynamic Properties.xcd和Ge Debye Temperature.xcd。 所示图形如下：
没有非谐性的实验结果表明在高温极限的Debye温度是395(3)K。 模拟的Debye温度是396 K，很好的与实验值相符。考虑到所完成的 计算的级别（given the level of calculation performed），这个结果还 是可以接受的。使用更好的k点格子可以提高精确度。通过实验， 推导出Debye温度在低温极限时的数值为374K，而CASTEP的预测 值为460K。在低温极限时的错误可通过下面的事实来解释，在G点 的频率并不是严格为0。 我们也应该注意到类似的谐波近似是声子计算的基础，但是当 温度高于Debye温度的三分之一时，类似的谐波近似是无效的。结 果，在温度高于Debye温度时，非谐性影响出现，不能够仅仅依靠 量的改变来解释。 尽管如此，总的来说，实验图表和由CASTEP产生的图表在质量 上是非常相似的。当Debye 最低温度达到255K时，在实验图表上会 骤然降低25K; CASTEP 显示出温度骤降出现在相同的位置并且预测 出Debye温度的最小值大约为270K。END