中科大-Materials-Studio-弹性常数计算.ppt

背景背景:当当前前，可可应应用用于于大大周周期期性性体体系系的的密密度度泛泛函函理理论论（DFTDFT）取取得得了了显显著著的的进进展展，已已经经成成为为解解决决材材料料设设计计、加加工工中中难难题题的的有有效效方方法法。人人们们依依据据这这个个理理论论可可以以使使解解释释实实验验数数据据，预预测测新新晶晶体体的的结结构构、结结合合能能和和表表面面活活性性等等基基本本性性质质。这这些些工工具具可可以以用用来来指指导导设设计计新新材材料料，允允许许研研究究人人员员理理解解基基本本的的化化学学和和物理过程。物理过程。绪论绪论:在在本本教教程程中中，将将学学习习如如何何使使用用CASTEPCASTEP来来计计算算弹弹性性常常数数和和其其他他的的力力学学性能。首先我们要优化性能。首先我们要优化BNBN立方晶体的结构，然后计算它的弹性常数。立方晶体的结构，然后计算它的弹性常数。本指南主要包括以下内容：本指南主要包括以下内容：1 1 优化优化BNBN立方晶体的结构立方晶体的结构 2 2 计算计算BNBN的弹性常数的弹性常数 3 3 弹性常数文件的描述弹性常数文件的描述使用使用CASTEPCASTEP计算计算BNBN的弹性常数的弹性常数目的：目的：使用使用 CASTEP CASTEP 计算弹性常数计算弹性常数模块：模块：Materials Visualizer,CASTEPMaterials Visualizer,CASTEP前提：前提：已使用已使用first principlesfirst principles预测了预测了AlAsAlAs的晶格常数的晶格常数1.1.优化优化BNBN立方晶体的结构立方晶体的结构 在计算弹性常数之前并不一定要进行几何优化，可以由实验观测到的结在计算弹性常数之前并不一定要进行几何优化，可以由实验观测到的结构计算出构计算出C Cijij数据。尽管如此，如果我们完成晶胞的几何优化，可以获得更多数据。尽管如此，如果我们完成晶胞的几何优化，可以获得更多相容的结果，进而计算与理论基态对应的弹性常数。相容的结果，进而计算与理论基态对应的弹性常数。弹性常数的精确度，尤其是切变常数的精确度，主要取决于弹性常数的精确度，尤其是切变常数的精确度，主要取决于SCFSCF计算的计算的品质，特别是布里渊区取样和波函数收敛程度的品质。所以我们设置品质，特别是布里渊区取样和波函数收敛程度的品质。所以我们设置SCFSCF、k k点取样和点取样和FFTFFT格子的精度为格子的精度为FineFine。首先导入首先导入BNBN结构结构 在在 菜菜 单单 栏栏 中中 选选 择择 File/File/ImportImport，从从structures/semiconductorsstructures/semiconductors中中选选中中BN.msiBN.msi，按按ImportImport按按钮钮，输输入入BNBN的晶体结构，见右图。的晶体结构，见右图。为了节省计算时间，由为了节省计算时间，由Build Build/Symmetry/Primitive Cell/Symmetry/Primitive Cell将此将此conventional conventional representation representation 转化为转化为primitive representation.primitive representation.现在设置几何优化现在设置几何优化 从工具栏中选择从工具栏中选择CASTEPCASTEP工具工具 ，然后从下拉列中选择然后从下拉列中选择CalculationCalculation（或从（或从菜单栏中选择菜单栏中选择Modules/CASTEP/Modules/CASTEP/CalculationCalculation）。）。CASTEP CalculationCASTEP Calculation对话框见右图：对话框见右图：在在SetupSetup标签中，把标签中，把TaskTask设置为设置为Geometry OptimizationGeometry Optimization，把，把Quality Quality 设设置为置为FineFine，并且把，并且把FunctionalFunctional设置为设置为GGA GGA and PW91and PW91。按下按下moremore按钮，选中按钮，选中Optimize cellOptimize cell。关闭关闭CASTEP Geometry OptimizationCASTEP Geometry Optimization对话对话框。框。选选择择ElectronicElectronic标标签签，按按下下More.More.按按钮钮以以得得到到CASTEP CASTEP Electronic Electronic OptionsOptions对对话话框框。把把Derived Derived gridgrid的的设设置置从从StandardStandard改改为为FineFine。关关闭闭CASTEP CASTEP Electronic OptionsElectronic Options对话框。对话框。选择选择Job ControlJob Control标签，设定本地机运算。标签，设定本地机运算。按下按下CASTEP CalculationCASTEP Calculation对话框中的对话框中的RunRun按钮。按钮。优化之后，此结构的晶胞参数应为优化之后，此结构的晶胞参数应为a=b=c=2.574a=b=c=2.574。现在我们可以继。现在我们可以继续计算优化结构的弹性常数。续计算优化结构的弹性常数。或按右键显示或按右键显示2.2.计算计算BNBN的弹性常数的弹性常数 BN CASTEP GeomOpt/BN.xsdBN CASTEP GeomOpt/BN.xsd处处于激活状态。于激活状态。选择选择CASTEP CASTEP CalculationCalculation对话框中的对话框中的SetupSetup标签，标签，从从TaskTask的下拉清单中选择的下拉清单中选择Elastic Elastic ConstantsConstants。按下按下More.More.按钮，按钮，CASTEP CASTEP Elastic ConstantsElastic Constants对话框见右图。对话框见右图。将将Number of steps for each Number of steps for each strainstrain由由4 4增加为增加为6 6，按，按RunRun运行。运行。CASTEPCASTEP的的弹弹性性常常数数计计算算任任务务的的结结果果以以一一批批.castep.castep输输出出文文件件的的形形式式给给出出。这这些些文文件件中中的的每每一一个个文文件件都都代代表表确确定定的的晶晶胞胞在在假假设设的的应应变变模模式式和和应应变变振振幅幅下下的的几几何何优优化化运运行行结结果果。这这些些文文件件的的命命名名约约定定为为：seedname_cij_m_nseedname_cij_m_n。对对于于给给定定的的模模式式来来说说，m m代代表表当当前前的的应应变变模模式式，n n代代表当前的应变振幅。表当前的应变振幅。6仅取一种应仅取一种应变模式变模式 从属性清单中选择从属性清单中选择Elastic constantsElastic constants，从，从BNBN的弹性常数计算工作中得到的结果文件的弹性常数计算工作中得到的结果文件BN.castepBN.castep应应自动显示在自动显示在Results fileResults file选框中。按下选框中。按下CalculateCalculate按按钮。计算结束后产生一个新的文档钮。计算结束后产生一个新的文档BN Elastic BN Elastic Constants.txtConstants.txt。此文档中的信息包括此文档中的信息包括：*输入的应变和计算出的应力的总结输入的应变和计算出的应力的总结*每一种应变模式线性拟合和拟合质量的计算结果每一种应变模式线性拟合和拟合质量的计算结果*给定对称性下计算出的应力与弹性常数之间的对应给定对称性下计算出的应力与弹性常数之间的对应*弹性常数弹性常数CijCij和弹性柔量和弹性柔量SijSij的表格的表格*导出量：体积模量和其倒数、压缩系数、杨氏模量、导出量：体积模量和其倒数、压缩系数、杨氏模量、PoissonPoisson比、比、Lame Lame 常数常数(用于模拟各向同性介质用于模拟各向同性介质)CASTEPCASTEP可可以以使使用用这这些些结结果果来来分分析析每每一一个个运运行行计计算算出出来来的的压压力力张张量量，产产生生一个有关弹性性质的文件。一个有关弹性性质的文件。从从工工具具栏栏中中选选择择CASTEP CASTEP 工工具具，然然后后选选择择AnalysisAnalysis或或者者从从菜菜单单栏栏中选择中选择Modules|CASTEP|AnalysisModules|CASTEP|Analysis。3 3 弹性常数文件的描述弹性常数文件的描述 对于这种点阵类型，需要考虑两种应变模式对于这种点阵类型，需要考虑两种应变模式(本教程只计算了一种本教程只计算了一种)。对于。对于每一种应变模式，都有一个计算出的应力的总结每一种应变模式，都有一个计算出的应力的总结(由各自的由各自的.castep.castep文件得到文件得到)。=Elastic constants from Materials Studio:CASTEP =Summary of the calculated stresses *Strain pattern:1=Current amplitude:1 Transformed stress tensor(GPa):-4.990578 0.000000 0.000000 0.000000 -6.907159 0.953658 0.000000 0.953658 -6.908215 Current amplitude:2 Transformed stress tensor(GPa):-5.949042 0.000000 0.000000 0.000000 -7.093625 0.571307 0.000000 0.571307 -7.094263 提供了应力，应变的组成和弹性常数张量之间联系的所有信息。在这一阶段，每一个弹性常数均有一个简洁的指数代表而不是由一对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 -4.990578 -0.003000 1 1 -5.949042 -0.001800 1 1 -6.891618 -0.000600 1 1 -7.838597 0.000600 1 1 -8.784959 0.001800 1 1 -9.726562 0.003000 C(gradient):788.920238 Error on C :0.945626 Correlation coeff:0.999997 Stress intercept:-7.363559 此梯度提供了弹性常数的数值（或弹性常数的线性组合），适配的质量，由相关系数表示，提供了另人满意的弹性常数的不确定度。在进一步的分析中没有使用压力的切点值,它很简单的指示出收敛的基态离最初的结构有多远。所有应变模式的结果总结如下：=Summary of elastic constants =id i j Cij(GPa)1 1 1 788.92024+/-0.946 4 4 4 447.55108+/-0.749 7 1 2 148.70983+/-0.754 The 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)=Elastic Compliance Constants Sij(1/GPa)=0.0014282 -0.0002075 -0.0002075 0.0000000 0.0000000 0.0000000 -0.0002075 0.0014282 -0.0002075 0.0000000 0.0000000 0.0000000 -0.0002075 -0.0002075 0.0014282 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0023533 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0023533 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0023533文件的最后部分包含推出的属性:736.57379 125.20883 125.20883 0.00000 0.00000 0.00000 125.20883 736.57379 125.20883 0.00000 0.00000 0.00000 125.20883 125.20883 736.57379 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 424.93974 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 424.93974 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 424.93974Bulk modulus =362.11330+/-0.593(GPa)Compressibility=0.00276(1/GPa)Axis Young Modulus Poisson Ratios (GPa)X 741.74894 Exy=0.1586 Exz=0.1586 Y 741.74894 Eyx=0.1586 Eyz=0.1586 Z 741.74894 Ezx=0.1586 Ezy=0.1586 Lame constants for isotropic material(GPa)Lambda=-106.1819,Mu=447.5511 END预测锗的热力学属性背景 线性响应或密度功能混乱理论是点阵动力学从头开始计算中最受欢迎的方法之一，尽管如此，这种方法的应用已经扩充到对振动属性的研究。线性响应提供了一种分析方法用于计算给定混乱的二级派生的整体能量。可以计算出许多属性，主要依赖于混乱的种类。在离子位置的混乱可以引起动力矩阵和声子；在磁场中引起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，把k-point 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并且把q-vector 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.500000 0.250000 0.750000)0.076923 +Effective cut-off=252.3267 EV +-+1 142.467048 +2 154.267057 +3 226.736315 +4 228.208353 +5 285.654435 +6 292.950626 +-+q-pt=2(0.500000 0.300000 0.700000)0.038462 +Effective cut-off=246.7136 EV +-+每一个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