晶格正交化[讨论]

卡开发发

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

这个过程中，相当于把两个晶格基矢线性组合，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），因此这样的情况其实只要建立一个合适大小的超胞就能够解决问题。

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

zhou

如果没有理解错误，即使满足上面的条件只有晶胞参数也不能还原原来晶胞，比如下图
左边格子与右边底边和高相同，按照规则正交化后的晶胞参数是一样的。

只有矩形晶胞的晶胞参数和原来晶胞的原子直角坐标，如何得到完整晶体呢？按照矩形晶胞的正交基矢平移得到的是上面图中左边两个A点，而实际晶体中并没有左上的A点。若果平移得到黄色的晶胞，则需要额外平移信息。

bashan

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

设所要 ...

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

端培阳

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

tjuptz

看神仙打架

卡开发发

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

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

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

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

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

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

granvia

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

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

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

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

4. 编程实现很容易啊。我曾写过简单程序，寻找overlayer和substrate之间的最小commensurate超胞，而且 允许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 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的情况可以放宽这种严格性，仅从模型研究的角度我同意你的观点，如果你要强调这种不失一般性，能写个实用性强程序分享给大家最好不过。

granvia

卡开发发 发表于 2019-7-9 21:31
“u和v的长度要么简单的有理关系”这点真的不一定，严格说“cos(t)虽不是有理数”和“只要满足|v|/|u|*co ...

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

舍入问题根本就不是问题：当你的精度足够严格时，总能找到一个足够大的正交原胞（尽管可能很大），其与原来的结构是完全等价的（在你所预订的精度内）。

granvia

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楼。

卡开发发

granvia 发表于 2019-7-9 21:53
严格说的确如此。但是第一性原理的很多模拟都是建立在这些commensurate近似的基础上的。最常见的例子就是 ...

其实上面做这样的讨论本来目的是看看是否有可能把非正交格子适当处理然后给Multiwfn分析，所以才考虑等价问题，否则直接挖团簇要简单得多，也可以按照特定的格点自己去写分析程序（虽然现在懒得干这样的事情）。对于异质结之类的overlayer，对于模型研究确实可以放宽要求，而且要比这个问题复杂，因为：最终的结构不见得要求正交还有你提到的旋转以及偏移。

granvia

卡开发发 发表于 2019-7-9 21:31
“u和v的长度要么简单的有理关系”这点真的不一定，严格说“cos(t)虽不是有理数”和“只要满足|v|/|u|*co ...

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

shgpei

granvia 发表于 2019-7-9 21:09
还需要知道你的a/c的比值大小

老师，我的a/c=1.7005 (晶胞参数:a=16.811 Å,b=4.7281 Å,c=9.886 Å,β=101.950°,α=γ=90°)，请问算的思路应该怎样呢？

卡开发发

granvia 发表于 2019-7-9 20:34
在理论上说，cos(t)虽不是有理数，只要满足|v|/|u|*cos(t)是有理数即可。但这样的coincidence几乎不太可 ...

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

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