表界面物性测试

刚刚开始接触MS计算，对表面能计算比较有兴趣，现在有以下问题求助各位大师：
1. 表面能我看可以用morphology或者CASTEP模块计算，这两者计算结果差别在哪里？我试着用morphology算了一下，出来的晶面的Esurf都是负数，这个照理说不合理啊？并且也没有单位。
2. 我看大家都在用CASTEP模块计算，这个建模过程中有如下细节不是很了解：
a. 单胞首先需要进行geometry optimization，LDA或者GGA算完，晶胞参数大概和实验值差别在什么正负百分比范围内可以接受呢？
b. 在cleave完surface之后加真空层，经常出现原子在真空层上面的情况，我查了一下cuturl('http://emuch.net/bbs/viewthread.php?tid=2870689')这个帖子里面有涉及，我也把这个改成Original就回去了，但是这个真的只是周期性排列的原因，对计算结果不会有影响吗？看着总是很别扭。
c. 模型建立以后我把下面几层固定了，然后上面几层弛豫，那下表面不和真空层接触的部分是否需要加H让键位饱和呢？还有这个超胞的其他侧面呢？
d. 文献上说LDA算表面比较靠谱，好像也是那个什么ladder里面最低一级的算起来也便宜点，不知道这个选择是否有差别。
3. 收敛之后这个具体表面能的计算方法：
a. Esurf = (Eslab  nEbulk)/A, where Esurf is the surface energy, Eslab is the total energy per repeated slab supercell, Ebulk is the energy per unit cell in the bulk, n is the number of unit cells and A is the total surface area per repeated unit.
我的理解就是把计算收敛结果status里面的Energy除以这个面的层数就是Eslab；而 A 可以比较方便的从构建的超胞的参数计算（A=a*b）。现在问题来了，这个Ebulk还有n是不是这么考虑：
把这个计算收敛的超胞的分子数搞出来，然后和每个单胞的分子数的比值就是n对吧？
b. 如果上面的考虑正确的话，哪位知道有什么方便的方法把超胞的分子式快速准确的读出来，而不是我一个个的去数吗？
4. 在计算表明能建模过程中，把下面几层固定，上面几层弛豫收敛，但是这个geometry optimization之后的能量是整个超胞的能量，我的想法是下面固定层能量会低一些，这样的话Eslab如果下面固定层越多岂不是表面的考虑越不准确？还请各位多多指教!

1. 表面能我看可以用morphology或者CASTEP模块计算，这两者计算结果差别在哪里？我试着用morphology算了一下，出来的晶面的Esurf都是负数，这个照理说不合理啊？并且也没有单位。
答：前者是分子力学的方法；后者是量子力学的方法，精度更高。Esurf为负值合理，若为正值，该表面不稳定，在现实中不会存在。单位应该是kcal/mol，建议仔细看看MS计算结果。
2. 我看大家都在用CASTEP模块计算，这个建模过程中有如下细节不是很了解：
a. 单胞首先需要进行geometry optimization，LDA或者GGA算完，晶胞参数大概和实验值差别在什么正负百分比范围内可以接受呢？
答:量子力学计算结果与实验值误差应≤2%。
b. 在cleave完surface之后加真空层，经常出现原子在真空层上面的情况，我查了一下cuturl('http://emuch.net/bbs/viewthread.php?tid=2870689')这个帖子里面有涉及，我也把这个改成Original就回去了，但是这个真的只是周期性排列的原因，对计算结果不会有影响吗？看着总是很别扭。
答：不影响结果。这完全是显示的问题，如果你能理解周期性边界条件，这个现象就很好理解。
c. 模型建立以后我把下面几层固定了，然后上面几层弛豫，那下表面不和真空层接触的部分是否需要加H让键位饱和呢？还有这个超胞的其他侧面呢？
答：固定的和能弛豫的晶面层，都属于晶体的结构，加H不是破坏了原结构么？个人认为不加，如果你研究的体系特殊，请参照相关文献来操作。
d. 文献上说LDA算表面比较靠谱，好像也是那个什么ladder里面最低一级的算起来也便宜点，不知道这个选择是否有差别。
答：通常LDA的精度没有GGA高，对于这些量子力学的方法，你可以都做计算，然后比较。
3. 收敛之后这个具体表面能的计算方法：
a. Esurf = (Eslab  nEbulk)/A, where Esurf is the surface energy, Eslab is the total energy per repeated slab supercell, Ebulk is the energy per unit cell in the bulk, n is the number of unit cells and A is the total surface area per repeated unit.
我的理解就是把计算收敛结果status里面的Energy除以这个面的层数就是Eslab；而 A 可以比较方便的从构建的超胞的参数计算（A=a*b）。现在问题来了，这个Ebulk还有n是不是这么考虑：
把这个计算收敛的超胞的分子数搞出来，然后和每个单胞的分子数的比值就是n对吧？
答：精确起见，不是读status里的数值，因为这里面的数值存在更新滞后的问题，最好的做法是几何优化后结果，计算但电能。n的理解应该没错。
b. 如果上面的考虑正确的话，哪位知道有什么方便的方法把超胞的分子式快速准确的读出来，而不是我一个个的去数吗？
答：在properties中，下拉有个各选项，其中有显示当前结构化学式的项目。
4. 在计算表明能建模过程中，把下面几层固定，上面几层弛豫收敛，但是这个geometry optimization之后的能量是整个超胞的能量，我的想法是下面固定层能量会低一些，这样的话Eslab如果下面固定层越多岂不是表面的考虑越不准确？
答：作收敛性测试，一般研究的晶面层数不能过少，但过多计算缓慢，作收敛性测试以确定既能保证精度又在合理的计算时间之内的层数。

帮主就是牛气！太感谢了，还有几点小问题没搞明白，再请教一下：
3. 收敛之后这个具体表面能的计算方法：
a. Esurf = (Eslab  nEbulk)/A, where Esurf is the surface energy, Eslab is the total energy per repeated slab supercell, Ebulk is the energy per unit cell in the bulk, n is the number of unit cells and A is the total surface area per repeated unit.
我的理解就是把计算收敛结果status里面的Energy除以这个面的层数就是Eslab；而 A 可以比较方便的从构建的超胞的参数计算（A=a*b）。现在问题来了，这个Ebulk还有n是不是这么考虑：
把这个计算收敛的超胞的分子数搞出来，然后和每个单胞的分子数的比值就是n对吧？
答：精确起见，不是读status里的数值，因为这里面的数值存在更新滞后的问题，最好的做法是几何优化后结果，计算但电能。
再问：第二次geometry optimization我没有选择optimize cell，status里面的energy值是收敛后的最终数据，这个也会有更新滞后的问题吗？更重要的是怎么计算单电能？我在.castep文件里面找了老半天也没有看到任何其他和能量相关的参数，也就是这个原因采用的status里面那个数值。是不是在calculation-property参数设置上有问题？还请明示。
b. 如果上面的考虑正确的话，哪位知道有什么方便的方法把超胞的分子式快速准确的读出来，而不是我一个个的去数吗？
答：在properties中，下拉有个各选项，其中有显示当前结构化学式的项目。
再问：这个也是同样的情况，我眼睛都看疼了也米找到这个玩意，也请再示意一下。对了，我用的是MS7.0。谢谢了！

再问：第二次geometry optimization我没有选择optimize cell，status里面的energy值是收敛后的最终数据，这个也会有更新滞后的问题吗？更重要的是怎么计算单电能？我在.castep文件里面找了老半天也没有看到任何其他和能量相关的参数，也就是这个原因采用的status里面那个数值。是不是在calculation-property参数设置上有问题？还请明示。
答：status中的能量数值有时候会更新滞后，所以为了保险起见，单独计算单点能；CASTEP中的单点能计算，就是把task选为Energy，对几何优化后的结果作计算。
再问：这个也是同样的情况，我眼睛都看疼了也米找到这个玩意，也请再示意一下。对了，我用的是MS7.0。谢谢了！
答：见附图。

补充两点看法：
关于2a和2d：以我经验。。。。LDA和GGA就能差50%以上。表面能的实验少之又少，另外，实验的精确度有多少，谁知道？我也记得哪里提了句LDA算表面能好，但是是不是有理有据，还有待考证。总的来说，LDA的算出来比GGA大不少，但是我看过的关于表面能的高档次文章，都是用GGA算的。
关于第4点：如果你截一个slab，slab上下两个surfaces是等价的，那么你完全可以optimize整个slab，得到的能量差除以面积2倍既是表面能。但是对于某些表面，特别是某些极化表面，你截出来的表面不能可能上下表面等价(一般一个是阳离子terminated，一个是阴离子terminated)，这样的话，你除以2倍面积，只能得出这两个不同极化表面的平均值，不能区分两个不同极化面之间的区别。这个时候就引入了你说的方法来解决这个问题。其实本身也就表面几层的原子会和体材料有差别，层数多了就趋近体材料了，所以只要层数够多，就不会不准确。

多谢补充关于LDA和DDA的问题，我也看了DDA算起来更加昂贵，一分钱一分货嘛，哈哈。
这个关于表面能的计算过程中的截slab，我研究过的体系里面的确有极性面的情况，也就是你说的上下表面一正一副的情况。我现在认为也是只能取个平均值了，因为如果强制把这个玩意都弄成上下表面都为同一类型电荷原子层的话，单胞的Stoichiometry就会发生变化，这个引入的差别就更加不好说了。
归根到底还是算起来太慢限制太多了，很多时候只能在合理的范围内尽量的近似好了。

GGA不是DDA ：）
极性表面的问题，感觉你理解的有点问题：肯定不能把上下表面都弄成同一类原子层，而是像你4里面说的，把slab底下几层固定住，只弛豫上面的几层。这样的话，固定住的部分看作是体材料，引入它主要是要考虑内部体材料对表面层原子的影响，而表面层的原子进行弛豫。这样的话，计算表面能的时候，需要能量差除以一倍的截面面积即可，因为体系只有一个表面了。如此也可以考虑不同极化面的区别了。你的slab的厚度当然要选的足够厚，固定层厚度和弛豫层厚度大约一半一半吧。具体的你可以搜搜关于ZnO极化表面计算的文章，是可以区分不同极化面的。

GGA，呵呵，已经被搞晕了：）
你说的这个很有意思，倒是要在这里再请教一下。
我现在对于计算时建的这个模型是这么理解的，在slab上（默认c方向）加上真空层之后，这个周期性的结构也就是在c方向上下表面之间引入真空层让原来的周期性在c轴方向断了好计算表面能，而a,b方向其实还是连续的所以也就是在最后考虑表面能和体相能能差的时候不需要考虑其他4个侧面面积的原因。
我后来又想了一下，这种情况下可以多建几个slab，固定中间几层不动，让上下几个表面同时弛豫，然后引入厚一点的真空层，这样就可以既考虑到极限面的上下正负原子差别，又不会改变原来的Stoichiometry。这样建模计算出来的表面积似乎更加合理。
而你提出的这个方法，我觉得虽然只是弛豫了一个方向上的原子，但是未弛豫另外一个c方向表面虽然是固定的，但是真空层的引入其实也让其成为了表面原子（前提条件是在这里没有另外加H），所以我觉得也还是需要把能量差除以两倍截面面积。
以上只是探讨，不足之处还请多多指教。

首先说你的多个slab的想法：我觉得是不成立的。对于极化面，你建N个slab，你就创造了N个阳极面和N个阴极面，所以还是看不出差别来。另外你要知道，特别对于极化面来说，slab要建的很厚的，一个面富集正电荷一个面富集负电荷，你想消除两个面之间的相互影响，需要的板和真空层厚度要比非极化面厚不少。
对于我说的方法：比如建一个slab，在中间（垂直于c轴的面）把板分成两个区域，一个叫surface（原子弛豫），一个叫bulk（原子固定）。surface即我们关心的表面。而引入bulk，只是为了取得内部体材料原子对表面原子的作用，所以两个区域都需要足够的厚度。弛豫过后，计算表面能的时候，体材料部分的能量是不计入的。surface的上表面暴露在真空中，是我们关心的，而下表面利用传统方法也是暴露在真空中的，但是用这个方法是与bulk链接的，即不为表面了。这个方法的缺点是，需要引入比常规方法多近一倍的原子数。对于极化面，刚才说了，本身板和真空就需要比较厚，所以算起来更贵了。

我看懂你的意思了，这个库伦作用力在上下表面两个极性面之间需要消除肯定会比非极性面需要更厚的真空层。
我现在唯一一点和你看法不同的就是这个下表面固定的bulk原子，你认为这个面是和bulk连接的，我认同其沿着c轴往上是周期性和bulk有bond,但是沿着c轴往下其实还是一个周期性的和上表面的真空层，所以我还是认为其是表面，只是我们的考量差别在于建模的时候是否让这个下表面弛豫。所以这也是在计算表面能的时候大家都是除以2A（A=a*b，默认沿着c轴弛豫），而不是只除以上表面的A。至于为什么惯例是只弛豫一个表面比如沿着c轴的上表面，我的看法是在大多数情况下如果不是极性面，上下表面差别不大，弛豫和已优化过的单胞固定的表面差别引入误差不会那么大。
在淡江大学李老师的slides里面找了个si的111真空层模型，我稍微改了一下，看看能不能说明这个问题。一个是单胞周期性结构，另外一个是我沿着c轴延展了一下，这样看起来下表面虽然在单胞里面看起来没有直接和真空层接触，但是正是因为这个玩意是周期性结构，其实一延展就能看见下表面其实是和其下方另一个单胞的上表面隔了一个真空层，这也是我理解的为什么要引入真空层以让周期性在c轴打破可以有表面弛豫的意义，以及需要改变其真空层厚度做收敛测试的原因，还请多指教。

Untitled-1.jpg

1.png