计算化学公社

标题: 晶格正交化[讨论] [打印本页]
作者
Author: 卡开发发    时间: 2017-4-6 06:10
标题: 晶格正交化[讨论]
本帖最后由 卡开发发 于 2017-4-7 13:47 编辑

特定情况下，我们希望晶体的晶格是正交的，如下图(图1)给出了2D晶格作为简单的例子，从白线到红线那样的方式做一个正交化。
(, 下载次数 Times of downloads: 101)

这个过程中，相当于把两个晶格基矢线性组合，v'=v，u'=c1*u+c2*v=|u|sinγ·e_v⊥。而更方便求出u'的方案就是先求出面积，A=u·v=|u||v|sinγ，然后除以v的模长|v|，得到长度|u'|，方向则是与v垂直的e_v⊥。因为体系是周期性的，所以只需要保证两晶格基矢相对位置不变即可，晶格基矢可以任意旋转平移。若推广到3D情形，第三个轴向w也类似，可以通过体积V=|u||v||w|√[1-cos^2(α)-cos^2(β)-cos^2(γ)+2cosαcosβcosγ]除以底面积求得模长，矢量则垂直于u'v'构成的平面；围成底面积的两基矢的正交化过程就和前面2D情形完全相同；原子的直角坐标还是不变的。
但上述的做法会出现9楼http://bbs.keinsci.com/forum.php ... d=42616&fromuid=308提到的问题，也就是说晶格内部的原子不能满足平移性，为了避免这样的问题，实际上我们可以再做进一步的处理：虽然按照这样的方式平移会在一个分量上偏离原有坐标，但对于特定的晶格偏离若干次就会回到原处（图2），因此这样的情况其实只要建立一个合适大小的超胞就能够解决问题。
(, 下载次数 Times of downloads: 101)

实际情形复杂很多，大部分的非正交晶格如果这样做正交化偏移的部分与相应的晶格参数之比根本就不是有理数，那么即便建立无穷大的超胞也不可能找到等价的正交晶格。此外，有些晶格正交的过程并不能限制v=v'，而是u'=c11u+c12v，v'=c21u+c22v，有些本应该能够正交化的晶格也就被排除了，因此原来的C++程序可能很多情况得不到我们要的结果，所以也只好撤回这个程序。对于二维情形，如果按照u'=c11u+c12v，v'=c21u+c22v考虑，以及考虑晶格内任意一个原子平移m*n次（两个维度分别平移）能够与原来重合，问题就变得颇为复杂，笔者暂时没想到有效的算法。虽然有一些小程序也能做，但反而不及手工做redefine lattice来的灵活，大家有啥好的想法欢迎讨论。



作者
Author: ruanyang    时间: 2017-4-6 13:47
抢个沙发，顶！
作者
Author: zhou    时间: 2017-4-6 20:30
这种正交化和直接写个基组（1,0,0），（0,1,0），（0，0，V）有什么区别吗
作者
Author: 卡开发发    时间: 2017-4-6 20:55
zhou 发表于 2017-4-6 20:30
这种正交化和直接写个基组（1,0,0），（0,1,0），（0，0，V）有什么区别吗

也可以直接那样写(u' 0 0) (0 v' 0) (0 0 w')，不过你得计算u' v' w'，并不是任意的长度都能满足条件。如果有更好的解决方案也欢迎讨论。
作者
Author: zhou    时间: 2017-4-6 21:07
卡开发发 发表于 2017-4-6 20:55
也可以直接那样写(u' 0 0) (0 v' 0) (0 0 w')，不过你得计算u' v' w'，并不是任意的长度都能满足条件。如 ...

如果只把晶胞参数改成正交的，加上原子的直角坐标并不能给出晶体正确的结构信息，还需要正交晶胞的周期性平移的矢量，这些矢量不是正交的。我在看gromacs手册的时候他也提到把晶胞变成砖块，这样做有什么实际的好处吗？
作者
Author: ggdh    时间: 2017-4-6 21:12
是不是任意的二维胞 都可以转换为正交胞啊
这样的话。做不同胞的复合层就简单多了。
作者
Author: 卡开发发    时间: 2017-4-6 21:16
zhou 发表于 2017-4-6 21:07
如果只把晶胞参数改成正交的，加上原子的直角坐标并不能给出晶体正确的结构信息，还需要正交晶胞的周期性 ...

肯定不是简单改成正交的，实际上严格一些的条件应该是要u'和u重合，v'在uv平面，实际这样做起来不复杂，这里只是提供一个思路。我上面提到了，以二维为例，如果要进行正交化又得保证给出晶体信息（当然对称性肯定是丢了）必须v=v'，u'=c1*u+c2*v=usinγ。

这样做没啥特别的目的，一来是近期有个朋友问我cp2k的电势格点拿到multiwfn怎么分析，我看看问题不算复杂，就打算帮忙解决一下而已；二来就是如BigDFT貌似至今只支持正交格子的计算，也许能帮上些忙，
作者
Author: 卡开发发    时间: 2017-4-6 21:36
ggdh 发表于 2017-4-6 21:12
是不是任意的二维胞 都可以转换为正交胞啊
这样的话。做不同胞的复合层就简单多了。

应该不是，比如二维上面的石墨烯的单胞也是是primitive cell，这种情况只能构造超胞再转换。还有没有其他特例我还没想出来。

目前来看很多做周期系统的人也急需Multiwfn的支持，总之欢迎大家讨论，能够提供更好的方案。
作者
Author: zhou    时间: 2017-4-6 21:41
如果没有理解错误，即使满足上面的条件只有晶胞参数也不能还原原来晶胞，比如下图 (, 下载次数 Times of downloads: 106)
左边格子与右边底边和高相同，按照规则正交化后的晶胞参数是一样的。

只有矩形晶胞的晶胞参数和原来晶胞的原子直角坐标，如何得到完整晶体呢？按照矩形晶胞的正交基矢平移得到的是上面图中左边两个A点，而实际晶体中并没有左上的A点。若果平移得到黄色的晶胞，则需要额外平移信息。
作者
Author: 卡开发发    时间: 2017-4-6 22:22
zhou 发表于 2017-4-6 21:41
如果没有理解错误，即使满足上面的条件只有晶胞参数也不能还原原来晶胞，比如下图
左边格子与右边底边和高 ...

有道理，我研究研究，明天搞个新的版本放上来。
作者
Author: Warm_Cloud    时间: 2017-4-7 10:10
本帖最后由 Warm_Cloud 于 2017-4-7 11:37 编辑

厉害
作者
Author: ggdh    时间: 2017-4-7 15:41
卡开发发 发表于 2017-4-6 21:36
应该不是，比如二维上面的石墨烯的单胞也是是primitive cell，这种情况只能构造超胞再转换。还有没有其他 ...

如果能构建出超胞再转，那也是很不错的啊
作者
Author: 卡开发发    时间: 2017-4-7 15:56
ggdh 发表于 2017-4-7 15:41
如果能构建出超胞再转，那也是很不错的啊

看样子恐怕是不行，9楼的说法是对的，不能只是做简单的平移。我把帖子改了，目前来说我没有想到有很好的算法。
作者
Author: 风之子    时间: 2017-9-15 12:11
一般我在处理的时候都是做切面，定义uv的值，通常六方的我会定义U（100）V（120），这样就得到正交的格子。对于其他角度的格子，会根据矢量的关系定义，一般都能得到正交化的格子。
作者
Author: 风之子    时间: 2017-9-15 12:14
风之子 发表于 2017-9-15 12:11
一般我在处理的时候都是做切面，定义uv的值，通常六方的我会定义U（100）V（120），这样就得到正交的格子。 ...

补充一下，正交后的一般会做超胞检查
作者
Author: 卡开发发    时间: 2017-9-15 16:32
风之子 发表于 2017-9-15 12:11
一般我在处理的时候都是做切面，定义uv的值，通常六方的我会定义U（100）V（120），这样就得到正交的格子。 ...

有想过这方面的问题，比如找最小公倍数之类的，六方格子正交化其实不是很困难，但确实会有可能会出现有些晶格正交化不了情况。看过其他资料，这个问题处理的也都不是很完美，所以只好换思路，比如通过FFT来求导之类的，能够做一些简单的电子密度方面的分析，暂时无其他进展。
作者
Author: xiaobudong    时间: 2017-9-21 20:50
风之子 发表于 2017-9-15 12:11
一般我在处理的时候都是做切面，定义uv的值，通常六方的我会定义U（100）V（120），这样就得到正交的格子。 ...

非六方的做不到绝对正交吧。
作者
Author: 风之子    时间: 2017-10-27 18:59
xiaobudong 发表于 2017-9-21 20:50
非六方的做不到绝对正交吧。

恩恩，不同的晶系有不同的取值
作者
Author: 风之子    时间: 2017-10-27 19:01
卡开发发 发表于 2017-9-15 16:32
有想过这方面的问题，比如找最小公倍数之类的，六方格子正交化其实不是很困难，但确实会有可能会出现有些 ...

是的，看似简单的反而没多少人研究得让人信服
作者
Author: shgpei    时间: 2019-7-9 15:36
风之子 发表于 2017-9-15 12:11
一般我在处理的时候都是做切面，定义uv的值，通常六方的我会定义U（100）V（120），这样就得到正交的格子。 ...

请问您可以具体说下根据矢量的关系定义，具体是怎么确定的呢？我的晶胞参数β=101.950°,α=γ=90°，望您不吝赐教
作者
Author: granvia    时间: 2019-7-9 19:11
本帖最后由 granvia 于 2019-7-9 22:34 编辑

可以简单地证明lz的想法是不可行的： 考虑一个2D原胞只有一个原子的情况，不妨将该原子放在原胞的一个顶点上（即原点，分数坐标(0,0)）。对于一般情况的平行四边形原胞，最邻近原子间的夹角不是直角。假设可以在不改变原胞面积的情况下将晶格矢量正交化，那么得到的正交原胞仍只含一个原子（因为原胞大小不变），但此时最邻近原子间的夹角却都是直角了，因此与之前的最邻近原子间的夹角不是直角的事实相矛盾。

当然上面考虑的是一般情况，一些例外是当原胞的一个基矢与两基矢之差正好正交，最简单的例子是原来的原胞平行四边形的内角是60度且两边长比值为1:2，此时最邻近原子间的夹角仍是直角。

构建超胞则有可以实现正交化，其算法其实很简单（以二维情况为例）：设超胞的晶格矢量为u'=m*u和v'=n*v-m*u。那么只需要找到u'和v'的夹角是90度即可：u'.v'=0，展开后可以导出如下关系：m/n = |v|*cos(t)/|u|，其中t是u与v的夹角。利用这个式子，我们先根据原来的原胞算出m/n的比值。显然，当且仅当cos(t)是有理数（根据Niven定理，t只能是60度或120度），m和n才有整数解。
求出的新基矢长度是：|u'| = m*|u|，|v'| = sqrt( n^2*|v|^2 - m^2*|u|^2 )


对于一般情况，即 cos(t)是无理数（即t是除60度和120度以外的任何角度时），那么我们可以利用辗转相除法（连分数法）得到最接近它的有理数，从而得到所需精度之内的最小超胞。

作者
Author: 卡开发发    时间: 2019-7-9 19:50
本帖最后由 卡开发发 于 2019-7-9 20:33 编辑
granvia 发表于 2019-7-9 19:11
可以简单地证明lz的想法是不可行的： 考虑一个2D原胞只有一个原子的情况，不妨将该原子放在原胞的一个顶点 ...

其实主要的问题是当时在做的时候光考虑了格子的正交，把里面原子的这件事忽略了，所以能够确定出发点肯定是错了，这点在前面几楼已经讨论得很清楚。后来也专门研究过了这方面的一些文献，有些算法确实可以做那样的近似正交化处理。不过总觉得即便能进行正交化也不是一个很讨巧的事情，所以也就没深入再研究了。不过，近似去把一个晶格做这样的变换（不一定非得正交）找最小超晶格还是有些用处，比如做异质结的时候。

还有就是，即便cos(t)不是有理数，其实u和v在特殊情况也能满足m/n是整数，反之，u和v在某些取值即便cos(t)是有理数正交化不见得有意义（实际可能需要很大很大的超晶格），u和本身也有舍入问题的，这样二维情况判断起来就可以很复杂，如果再考虑三维情形可能还要更复杂一些。

作者
Author: granvia    时间: 2019-7-9 20:34
本帖最后由 granvia 于 2019-7-9 20:40 编辑
卡开发发 发表于 2019-7-9 19:50
其实主要的问题是当时在做的时候光考虑了格子的正交，把里面原子的这件事忽略了，所以能够确定出发点肯定 ...

在理论上说，cos(t)虽不是有理数，只要满足|v|/|u|*cos(t)是有理数即可。但这样的coincidence几乎不太可能，因为u和v的长度要么简单的有理关系（如纯金属的晶体），要么都是不规则的（如分子晶体）。你说的“u和v在某些取值即便cos(t)是有理数也不见得能正交化”，我也不同意，因为只要给定精度，任何小数都可以用连分数法求出其有理数比值，你需要的精度即便再小，大不了我给的超胞足够大就行。如无理数pi可表达为 22/7， 333/106， 355/113，52163/16604，... ，可按照你所需精度要求，不断逼近下去即可。

没错，这样做的代价就是超胞可能非常非常大，这样就失去了实用意义。我这里只是给出这个问题的一般性算法，仅供各位参考。
作者
Author: granvia    时间: 2019-7-9 21:09
shgpei 发表于 2019-7-9 15:36
请问您可以具体说下根据矢量的关系定义，具体是怎么确定的呢？我的晶胞参数β=101.950°,α=γ=90°，望 ...

还需要知道你的a/c的比值大小
作者
Author: 卡开发发    时间: 2019-7-9 21:31
本帖最后由 卡开发发 于 2019-7-9 21:40 编辑
granvia 发表于 2019-7-9 20:34
在理论上说，cos(t)虽不是有理数，只要满足|v|/|u|*cos(t)是有理数即可。但这样的coincidence几乎不太可 ...

“u和v的长度要么简单的有理关系”这点真的不一定，严格说“cos(t)虽不是有理数”和“只要满足|v|/|u|*cos(t)是有理数这是两件事，只是因为数值舍入才导致一些无理数看起来像是有理数。舍入问题放宽标准其实意味着正交构造出来的晶体其实和原来的结构并不等价，否则前面也就没必要讨论那么多了。当然从前面来看，不见得所有的晶格结构都能找得到。

如果没有等价的限制，那总是能找到一些近似的正交结构，只是有可能超晶格非常大而已。

作者
Author: shgpei    时间: 2019-7-9 21:32
granvia 发表于 2019-7-9 21:09
还需要知道你的a/c的比值大小

老师，我的a/c=1.7005 (晶胞参数:a=16.811 Å,b=4.7281 Å,c=9.886 Å,β=101.950°,α=γ=90°)，请问算的思路应该怎样呢？
作者
Author: granvia    时间: 2019-7-9 21:53
卡开发发 发表于 2019-7-9 21:31
“u和v的长度要么简单的有理关系”这点真的不一定，严格说“cos(t)虽不是有理数”和“只要满足|v|/|u|*co ...

严格说的确如此。但是第一性原理的很多模拟都是建立在这些commensurate近似的基础上的。最常见的例子就是表面上overlayer体系的模拟，当两种表面的晶格周期性不匹配时（往往如此；而且当overlayer还有取向旋转时，还会出现moiré图样，则更加复杂了），只能在理论上找出精度允许范围内的最小公倍数超胞进行建模即可。这些都只是理论模型，理论模型虽不能100%还原现实，但能说明问题，揭示物理本质即可。
作者
Author: 卡开发发    时间: 2019-7-9 22:14
granvia 发表于 2019-7-9 21:53
严格说的确如此。但是第一性原理的很多模拟都是建立在这些commensurate近似的基础上的。最常见的例子就是 ...

其实上面做这样的讨论本来目的是看看是否有可能把非正交格子适当处理然后给Multiwfn分析，所以才考虑等价问题，否则直接挖团簇要简单得多，也可以按照特定的格点自己去写分析程序（虽然现在懒得干这样的事情）。对于异质结之类的overlayer，对于模型研究确实可以放宽要求，而且要比这个问题复杂，因为：最终的结构不见得要求正交还有你提到的旋转以及偏移。
作者
Author: granvia    时间: 2019-7-9 22:32
本帖最后由 granvia 于 2019-7-10 21:43 编辑
shgpei 发表于 2019-7-9 21:32
老师，我的a/c=1.7005 (晶胞参数:a=16.811 Å,b=4.7281 Å,c=9.886 Å,β=101.950°,α= ...

设新基矢量a'=m*a, c'=n*c-m*a，利用上述公式：m/n = |c|/|a|*cos(β) = 9.886/16.811*cos(101.950°) = -0.12176405...
这个比值可用分数表示： 1/8 = 0.125，5/41 = 0.12195...，14/115 = 0.12173...，33/271 = 0.121764...

由于你给的晶胞参数有效数字时4位，所以33/271足以精确表示你的超胞了（you can't do any better）。相应地，a' = 33*a，c' = 271*c + 33*a。这样，|a'| = 554.763，|c'| = 4544.1。 可见，这个超胞非常大。

如果采取简化模型，可取m/n = -1/8 ，则a' = a，c' = 8*c + a，|c'| = 134.1。

这个例子的更好结果及其算法和代码实现，见35楼。

作者
Author: granvia    时间: 2019-7-10 00:07
本帖最后由 granvia 于 2019-7-10 00:09 编辑
卡开发发 发表于 2019-7-9 21:31
“u和v的长度要么简单的有理关系”这点真的不一定，严格说“cos(t)虽不是有理数”和“只要满足|v|/|u|*co ...

“u和v的长度要么简单的有理关系”，最简单的就是你前面所举的石墨烯的例子，正是因为它的u和v比值是1:1，而且正好原胞的内角正好是60和120度，所以才能找到整数解的m和n。（但一般我们不把如石墨烯这样的六方晶胞进行正交化，否则就无法体现出其旋转对称性了）

舍入问题根本就不是问题：当你的精度足够严格时，总能找到一个足够大的正交原胞（尽管可能很大），其与原来的结构是完全等价的（在你所预订的精度内）。
作者
Author: 卡开发发    时间: 2019-7-10 07:59
本帖最后由 卡开发发 于 2019-7-10 08:20 编辑
granvia 发表于 2019-7-10 00:07
“u和v的长度要么简单的有理关系”，最简单的就是你前面所举的石墨烯的例子，正是因为它的u和v比值是1:1 ...

六方的格子|u|/|v|只是很特殊的例子，并没有普遍意义。另外，这种体系不进行正交化只是一般情况下没必要，与其旋转对称性无关，比如研究特定覆盖度的时候还是会进行正交化（比如Cu(1 1 1)、Al(1 1 1)的吸附等，这方面例子太多了）。

还有上面说到看cos(t)的说法其实不是很靠谱，比如你拿六方格子60度夹角的两个矢量做1x2的晶格，再选取u'=u，v'=u+2v，这时候夹角就是45°，你可以这样变换任意次，所得到的夹角t对应的cos(t)基本上应该都不是有理数，但这样的格子都能严格正交化。

涉及到舍入阈值的问题，正交化之后的格子也就无法说与原来的格子严格等价，尽管晶格内的原子化学环境确实与原来的结构接近，从实际应用角度来说，如果需要非常大的格子，这种情况对于我的需求而言不如使用团簇模型。这方面我觉得也没必要再争论下去。

关于构造overlayer的情况可以放宽这种严格性，仅从模型研究的角度我同意你的观点，如果你要强调这种不失一般性，能写个实用性强程序分享给大家最好不过。
作者
Author: SigFig    时间: 2019-7-10 09:06
正交，长方和正方类基本等于白给，六方类也还好处理。

对于单斜点群，其正交化的晶格面积应该比其最小晶格大。
最小面积什么样我不知道，如果有互质整数M和N使得M|a|=N|b|，那么可以得到一个面积为2MN倍于单晶胞面积的正交晶胞。其坐标为u=Ma+Nb,v=Ma-Nb (a,b为向量),夹角为90度。

如果a,b是无理数，这个晶胞可能会不存在。
在有限精度的情况下，一定存在，只是会很大。
注意一点，这个晶格的存在和夹角无关。

特殊情况，当点群为六方时，|a|=|b|，M=N=1，得到的就是以晶胞对角线形成的正交晶胞。

如果只是需要一个正交晶胞了话，3d体系也可以延续这个思路，只是对于三斜晶系晶胞会更大。
作者
Author: granvia    时间: 2019-7-10 16:24
卡开发发 发表于 2019-7-10 07:59
六方的格子|u|/|v|只是很特殊的例子，并没有普遍意义。另外，这种体系不进行正交化只是一般情况下没必要 ...

1. 六方格子之所以不转化为正交格子，是为了保留晶胞的对称性，六方的空间群对称元素显然要多于一般的正交格子——其中就包括旋转对称性。

2. 如果我对你的描述没有理解错的话，这里u'和v'的夹角不是45度，而是acos(2/sqrt(7))，约41度。

3. 这个我前面已经也说过了，并无异议。我这里在纯理论上进行探讨，仅供参考。

4. 编程实现很容易啊。我曾写过简单程序，寻找overlayer和substrate之间的最小commensurate超胞，而且 允许overlayer取向的旋转。
作者
Author: granvia    时间: 2019-7-10 16:33
本帖最后由 granvia 于 2019-7-10 16:52 编辑
SigFig 发表于 2019-7-10 09:06
正交，长方和正方类基本等于白给，六方类也还好处理。

对于单斜点群，其正交化的晶格面积应该比其最小晶 ...

你的这个算法很优雅，而且更具一般性。看来我上面的算法有很大局限性。 因此，M和N有无有理数解的确与a和b的夹角无关。

作者
Author: granvia    时间: 2019-7-10 17:19
本帖最后由 granvia 于 2019-7-10 21:41 编辑

参考了32楼和21楼的两种算法（它们之间并不等价或包含），我觉得可以提出更一般性的算法，如下：

设所要寻找的正交超胞的基矢为：u = m1*a + n1*b，v = m2*a + n2*b，则根据 u.v = 0 可得：m1*m2 + n1*n2*s + (m1*n2 + m2*n1)*s*cos(t) = 0，其中 s = |b|/|a|，t为a与b夹角。这个超胞的面积是原来晶胞的(m1*n2-m2*n1)倍。

当m1 = m2 = M，n1 = N，n2 = -N时，得到32楼提出的解；当m1 = m， m2 = -m，n1 = 0， n2 = n时，得到21楼提出的解。 但显然，存在更多可能的(m1,m2,n1,n2)组合，满足上述方程。 在给定精度下，这个四参数的近似解(m1,m2,n1,n2)有可能比32楼的(M,N)和21楼(m,n)解给出更小的超胞。换句话说，利用这个更一般性的方程，可以得到给定精度下的最优解。在算法上可以枚举搜索所有的m1，m2，n1，n2的整数组合，在满足精度要求下，找到使得(m1*n2-m2*n1)为最小正值的组合。至于精度，可以定义一个误差函数为：
  f = [ m1*m2 + n1*n2*s + (m1*n2 + m2*n1)*s*cos(t) ] / sqrt(m1^2 + n1^2*s^2 + 2*m1*n1*s*cos(t)) / sqrt( m2^2 + n2^2*s^2 + 2*m2*n2*s*cos(t) )
其中多出来的两个分母sqrt(...)是将u和v归一化的结果。据此，我们只需搜索所有满足 f <= err （精度）且使(m1*n2-m2*n1)为最小正值的(m1,m2,n1,n2)组合即可。

将这个算法应用于26楼的例子，结果如下（A表示超胞大小是原来晶胞的A倍）：
当精度为0.1时，m1 = 1; m2 = 1; n1 = -2; n2 = 1; A = 3; f = -3.1E-02
当精度为0.01时，m1 = 0; m2 = 5; n1 = -1; n2 = 1; A = 5; f = 7.2E-03
当精度为0.001时，m1 = 1; m2 = 5; n1 = -4; n2 = 3; A = 23; f = 9.7E-04


当精度为0.0001时，m1 = 3; m2 = -5; n1 = 2; n2 = 17; A = 61; f = 5.4E-05
（。。。我们可以将精度任意提高下去。。。）

可见，以上算法给出的最优超胞都不是32楼或21楼能给出的解，因此更具有一般性。当精度要求不高（0.01）时，只需扩胞5倍就能可以得到近似模型了。

具体MATLAB代码如下（为简便起见，只打印出一种可能的最优解）：
  1. function findOthroCell()

  3. % a=16.811 Angstrom, b=4.7281 Angstrom, c=9.886 Angstrom, beta=101.950, alpha=gamma=90 Deg.
  4. a = 16.811;
  5. b = 9.886;
  6. t = 101.950/180*pi;

  8. Err = 1e-4;
  9. Mmax = 10;
  10. Nmax = 10;

  12. s = b/a;
  13. st = s*cos(t);

  15. Amin = Inf;
  16. for m1 = 0 : Mmax
  17.     for n1 = -Nmax : Nmax
  18.         % a' = m1*a + n1*b != 0
  19.         if( m1 == 0 && n1 == 0 )
  20.             continue;
  21.         end
  22.         
  23.         % Length of new basis vector a:
  24.         len1 = sqrt( m1^2 + n1^2*s^2 + 2*m1*n1*st );
  25.         
  26.         for m2 = -Mmax : Mmax
  27.             for n2 = -Nmax : Nmax
  28.                 % b' = m2*a + n2*b != 0
  29.                 if( m2 == 0 && n2 == 0 )
  30.                     continue;
  31.                 end
  32.                
  33.                 % Length of new basis vector b:
  34.                 len2 = sqrt( m2^2 + n2^2*s^2 + 2*m2*n2*st );
  35.         
  36.                 % Precision:
  37.                 f = (m1*m2 + n1*n2*s + (m1*n2 + m2*n1)*st)/len1/len2;
  38.                 if( abs(f) > Err )
  39.                     continue
  40.                 end
  41.                
  42.                 A = m1*n2-m2*n1; % Area in units of original cell area
  43.                 if( A <= 0 )
  44.                     continue;
  45.                 end
  46.                 fprintf( 'm1 = %4i; m2 = %4i; n1 = %4i; n2 = %4i; ', ...
  47.                     m1, m2, n1, n2 );
  48.                 fprintf( ' err = %10.1E;  A = %12.6f  Amin = %12.6f', ...
  49.                     f, A, Amin );
  50.                 if( A < Amin )
  51.                     Amin = A;
  52.                     m1_opt = m1;
  53.                     n1_opt = n1;
  54.                     m2_opt = m2;
  55.                     n2_opt = n2;
  56.                     f_opt = f;
  57.                     fprintf( '(NEW Amin)' );
  58.                 end
  59.                 fprintf( '\n' );
  60.             end
  61.         end
  62.     end
  63. end

  65. fprintf( 'Optimal supercell within error of %.1E:\n', Err );
  66. fprintf( 'm1 = %i; m2 = %i; n1 = %i; n2 = %i; A = %i; f = %.1E\n', ...
  67.     m1_opt, m2_opt, n1_opt, n2_opt, Amin, f_opt );

  69. end
复制代码

作者
Author: 卡开发发    时间: 2019-7-11 00:59
granvia 发表于 2019-7-10 16:24
1. 六方格子之所以不转化为正交格子，是为了保留晶胞的对称性，六方的空间群对称元素显然要多于一般的正 ...

1、除非是原胞计算强调其对称性作用的情况下才会保留六方格子，而且不仅限于六方格子本身，其余情况不一定。
2、确实笔误，应该是41°，并不影响结果的讨论。
作者
Author: tjuptz    时间: 2019-7-11 14:25
看神仙打架
作者
Author: 端培阳    时间: 2023-9-28 18:13
您好！我想请教在用MS切表面时，怎么修改UV来使切出来的面小一点！诚恳求教！（抱歉我没权限发短消息，只能在这里留言！）（vx：wx2212509542）
作者
Author: bashan    时间: 2024-6-4 11:27
本帖最后由 bashan 于 2024-6-4 12:26 编辑
granvia 发表于 2019-7-10 17:19
参考了32楼和21楼的两种算法（它们之间并不等价或包含），我觉得可以提出更一般性的算法，如下：

设所要 ...

这个程序不错，如果把误差评价标准改成角度差，给出的结果将更加完美



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3