RASPA批量计算输入文件中的Unitcells值怎么确定 - 计算化学公社

import math
def get_unit_cell_number(framework_path, r_cut):
# 从cif文件中读取晶格常数
flag = [
'_cell_length_a', '_cell_length_b', '_cell_length_c',
'_cell_angle_alpha', '_cell_angle_beta', '_cell_angle_gamma'
]
flag_length = len(flag)
with open(framework_path, 'r') as f:
i = 1
flag_count = 0
while i <= 5e4:
flag_content = flag[flag_count]
line = f.readline()
if flag_content in line:
if flag_content == '_cell_length_a':
a = float(line.strip().split()[1])
elif flag_content == '_cell_length_b':
b = float(line.strip().split()[1])
elif flag_content == '_cell_length_c':
c = float(line.strip().split()[1])
elif flag_content == '_cell_angle_alpha':
alpha = float(line.strip().split()[1]) * math.pi / 180 # 计算时用弧度单位
elif flag_content == '_cell_angle_beta':
beta = float(line.strip().split()[1]) * math.pi / 180 # 计算时用弧度单位
elif flag_content == '_cell_angle_gamma':
gamma = float(line.strip().split()[1]) * math.pi / 180 # 计算时用弧度单位
else:
print(f'No specific operation for flag_content: {flag_content}!')
flag_count += 1
if flag_count == flag_length:
break
i += 1
# 计算晶胞体积
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
# 计算晶胞各个面的表面积
base_area_x = b * c * math.sin(alpha)
base_area_y = a * c * math.sin(beta)
base_area_z = a * b * math.sin(gamma)
# 计算各个方向的最小距离，即平行六面体各个方向的高，等于体积除以底面积
perpendicular_length_x = V / base_area_x
perpendicular_length_y = V / base_area_y
perpendicular_length_z = V / base_area_z
# 根据截断半径r_cut计算所需各方向的unit_cell数目
unit_cell_x = math.ceil(2 * r_cut / perpendicular_length_x)
unit_cell_y = math.ceil(2 * r_cut / perpendicular_length_y)
unit_cell_z = math.ceil(2 * r_cut / perpendicular_length_z)
return ' '.join([str(unit_cell_x), str(unit_cell_y), str(unit_cell_z)])

复制代码