qlx 发表于 2023-9-30 16:09
盒子需要大于两倍的cutoff，指的是盒子三个方向的两平行表面间的距离都需要大于两倍的cutoff。原代码计算的 ...

好的，非常感谢，确实可以解决问题。之前没有想到用体积进行换算

盒子需要大于两倍的cutoff，指的是盒子三个方向的两平行表面间的距离都需要大于两倍的cutoff。原代码计算的一条边向另外一条边的投影，并不能准确地反映非正交晶胞中两平行表面间的距离。以你附件中的CORE2019MOF484.cif为例，按原代码中的算法，会认为前后两个平行面间的距离是OA或AD中较小的那个数值（这里是AD），但实际上，无论是OA还是AD都不是OBC这个平面的法线方向，前后两个表面的距离，既不是OA也不是AD。



可以根据平行六面体的体积公式，由晶格常数算出体积，除以各个方向的底面积，得到高，即平行表面间的距离。平行六面体体积公式可参见维基百科。



示例代码如下（python）：

1. import math
   
2. 
   
3. def get_unit_cell_number(framework_path, r_cut):
   
4. 
   
5.     # 从cif文件中读取晶格常数
   
6.     flag = [
   
7.         '_cell_length_a', '_cell_length_b', '_cell_length_c',
   
8.         '_cell_angle_alpha', '_cell_angle_beta', '_cell_angle_gamma'
   
9.     ]
   
10.     flag_length = len(flag)
    
11.     with open(framework_path, 'r') as f:
    
12.         i = 1
    
13.         flag_count = 0
    
14.         while i <= 5e4:
    
15.             flag_content = flag[flag_count]
    
16.             line = f.readline()
    
17.             if flag_content in line:
    
18.                 if flag_content == '_cell_length_a':
    
19.                     a = float(line.strip().split()[1])
    
20.                 elif flag_content == '_cell_length_b':
    
21.                     b = float(line.strip().split()[1])
    
22.                 elif flag_content == '_cell_length_c':
    
23.                     c = float(line.strip().split()[1])
    
24.                 elif flag_content == '_cell_angle_alpha':
    
25.                     alpha = float(line.strip().split()[1]) * math.pi / 180 # 计算时用弧度单位   
    
26.                 elif flag_content == '_cell_angle_beta':
    
27.                     beta = float(line.strip().split()[1]) * math.pi / 180   # 计算时用弧度单位
    
28.                 elif flag_content == '_cell_angle_gamma':
    
29.                     gamma = float(line.strip().split()[1]) * math.pi / 180 # 计算时用弧度单位
    
30.                 else:
    
31.                     print(f'No specific operation for flag_content: {flag_content}!')
    
32.                
    
33.                 flag_count += 1
    
34.                 if flag_count == flag_length:
    
35.                     break
    
36.             i += 1
    
37. 
    
38.     # 计算晶胞体积         
    
39.     V =  a * b * c * (1 + 2 * math.cos(alpha) * math.cos(beta) * math.cos(gamma) - (math.cos(alpha))**2 - (math.cos(beta))**2 - (math.cos(gamma))**2) ** 0.5
    
40.    
    
41.     # 计算晶胞各个面的表面积
    
42.     base_area_x = b * c * math.sin(alpha)
    
43.     base_area_y = a * c * math.sin(beta)
    
44.     base_area_z = a * b * math.sin(gamma)
    
45. 
    
46.     # 计算各个方向的最小距离，即平行六面体各个方向的高，等于体积除以底面积
    
47.     perpendicular_length_x = V / base_area_x
    
48.     perpendicular_length_y = V / base_area_y
    
49.     perpendicular_length_z = V / base_area_z
    
50. 
    
51.     # 根据截断半径r_cut计算所需各方向的unit_cell数目
    
52.     unit_cell_x = math.ceil(2 * r_cut / perpendicular_length_x)
    
53.     unit_cell_y = math.ceil(2 * r_cut / perpendicular_length_y)
    
54.     unit_cell_z = math.ceil(2 * r_cut / perpendicular_length_z)
    
55. 
    
56.     return ' '.join([str(unit_cell_x), str(unit_cell_y), str(unit_cell_z)])

复制代码

slxc920113 发表于 2023-9-15 09:22
用atomsk批量对cif最正交化变换。

不是所有晶胞都可以正交变换的

用atomsk批量对cif最正交化变换。