晶格正交化[讨论]

卡开发发

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的情况可以放宽这种严格性，仅从模型研究的角度我同意你的观点，如果你要强调这种不失一般性，能写个实用性强程序分享给大家最好不过。

SigFig

正交，长方和正方类基本等于白给，六方类也还好处理。

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

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

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

如果只是需要一个正交晶胞了话，3d体系也可以延续这个思路，只是对于三斜晶系晶胞会更大。

granvia

卡开发发 发表于 2019-7-10 07:59
六方的格子|u|/|v|只是很特殊的例子，并没有普遍意义。另外，这种体系不进行正交化只是一般情况下没必要 ...

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

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

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

4. 编程实现很容易啊。我曾写过简单程序，寻找overlayer和substrate之间的最小commensurate超胞，而且 允许overlayer取向的旋转。

granvia

SigFig 发表于 2019-7-10 09:06
正交，长方和正方类基本等于白给，六方类也还好处理。

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

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

granvia

参考了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()
   
2. 
   
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;
   
7. 
   
8. Err = 1e-4;
   
9. Mmax = 10;
   
10. Nmax = 10;
    
11. 
    
12. s = b/a;
    
13. st = s*cos(t);
    
14. 
    
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
    
64. 
    
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 );
    
68. 
    
69. end

复制代码

卡开发发

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

1、除非是原胞计算强调其对称性作用的情况下才会保留六方格子，而且不仅限于六方格子本身，其余情况不一定。
2、确实笔误，应该是41°，并不影响结果的讨论。

tjuptz

看神仙打架

端培阳

您好！我想请教在用MS切表面时，怎么修改UV来使切出来的面小一点！诚恳求教！（抱歉我没权限发短消息，只能在这里留言！）（vx：wx2212509542）

bashan

granvia 发表于 2019-7-10 17:19
参考了32楼和21楼的两种算法（它们之间并不等价或包含），我觉得可以提出更一般性的算法，如下：

设所要 ...

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