From 78244319aa8addce24c194300f13d4895ecdb38c Mon Sep 17 00:00:00 2001 From: zqxian <2821838490@qq.com> Date: Sun, 22 Jun 2025 21:46:18 +0800 Subject: [PATCH 1/2] =?UTF-8?q?20250622=E4=BF=9D=E7=95=99=E5=90=B8?= =?UTF-8?q?=E9=99=84=E4=BD=8D=E7=82=B9=E5=92=8C=E5=90=B8=E9=99=84=E6=A8=A1?= =?UTF-8?q?=E5=BC=8F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- HTMACat/catkit/gen/adsorption.py | 93 +++++++++++++++++++++++--------- 1 file changed, 68 insertions(+), 25 deletions(-) diff --git a/HTMACat/catkit/gen/adsorption.py b/HTMACat/catkit/gen/adsorption.py index 14792e9..08fedaf 100644 --- a/HTMACat/catkit/gen/adsorption.py +++ b/HTMACat/catkit/gen/adsorption.py @@ -571,6 +571,7 @@ def get_principal_axes(inertia_tensor): sorted_indices = np.argsort(eigenvalues)[::-1] # 按惯性矩大小排序（从大到小） principal_axes = eigenvectors[:, sorted_indices] return principal_axes + def _single_adsorption( self, adsorbate, @@ -686,7 +687,66 @@ def _single_adsorption( # 根据参与吸附的原子确定位向，将物种“扶正” adsorption_vector, flag = utils.solve_normal_vector_linearsvc(atoms.get_positions(), bond) atoms.rotate(adsorption_vector, [0, 0, 1]) + elif direction_mode == 'asphericity': + # 根据分子的形状调整朝向，使其“平躺”在表面上 + masses = atoms.get_masses() + positions = atoms.get_positions() + #xzq + # 计算惯性矩张量 + inertia_tensor_value = self.inertia_tensor(positions, masses) + + # 获取主惯性轴 + principal_axes = self.get_principal_axes(inertia_tensor_value) + + # 自动处理三种惯性矩方向 + slabs_list = [] + for inertia_mode in range(1, 4): # 遍历第一、第二、第三惯性矩 + atoms_copy = atoms.copy() # ✅ 每次都从原始结构复制 + base_position = [0.0, 0.0, 0.0] + adsorption_vector = principal_axes[:, inertia_mode - 1] + atoms_copy.rotate(adsorption_vector, [0, 0, 1]) + + if enable_rotate_xoy and rotation_mode == 'vnn' and rotation_args != {}: + principle_axe = utils.solve_principle_axe_pca(atoms_copy.get_positions()) + if abs(rotation_args['vec_to_neigh_imgsite'][0]) < 1e-8: + target_vec = [1, 0, 0] + elif abs(rotation_args['vec_to_neigh_imgsite'][1]) < 1e-8: + target_vec = [0, 1, 0] + else: + target_vec = [-1/rotation_args['vec_to_neigh_imgsite'][0], + 1/rotation_args['vec_to_neigh_imgsite'][1], 0] + atoms_copy.rotate([principle_axe[0], principle_axe[1], 0], target_vec) + # 设置 吸附 高度 + z_coordinates = atoms_copy.get_positions()[:, 2] + min_z = np.min(z_coordinates) + + final_positions = slab.get_positions() + max_z = np.max(final_positions[:, 2]) + + + if abs(site_coord[2]) < 1e-6: # 如果为 0 或非常接近 0 + effective_distance = 2 + else: + effective_distance = site_coord[2] + + base_position[2] = round(max_z - min_z + effective_distance, 2) + + #使用输入的 site_coord 作为 xy 坐标 + base_position[0] = round(site_coord[0], 1) + base_position[1] = round(site_coord[1], 1) + atoms_copy.translate(base_position) # 确保所有构型都以正确的 z 轴间距平移 + + # 将分子平移到指定的 site_coord 位置 + vec_tmp = site_coord - atoms_copy.get_positions()[bond] # 计算平移量 + atoms_copy.translate([vec_tmp[0], vec_tmp[1], 0]) # 只平移 x-y，不改变 z + + + slab_with_adsorbate = slab.copy() + slab_with_adsorbate += atoms_copy + slabs_list.append(slab_with_adsorbate) + + return slabs_list # 返回三种吸附构型的列表) elif direction_mode == 'hetero': # zjwang 20240815 # 适用于有明确“官能团”的偏链状的分子 # 求出从分子杂原子中心到所有原子中心的吸附矢量vec_p，将vec_p旋转到[0,0,1]竖直向上 @@ -741,24 +801,7 @@ def _single_adsorption( a = copy.deepcopy(a_orig) # 确保每次处理原始未改动的 adsorbate dealt_positions = a.get_positions() min_z = np.min(dealt_positions[:,2]) #得到分子z轴最小值 - # Step 1: Z 方向平移 - z_translation = max_z - min_z + (2.5 if site_coord[2] == 0 else site_coord[2]) - a.translate([0, 0, z_translation]) - # 计算 slab 中心坐标 - slab_positions = slab.get_positions() - center_x, center_y = utils.center_slab(slab_positions) - # xzq如果输入的 xy 为 [0, 0]，则将分子放在 slab 的 xy 中心 - if site_coord[0] == 0 and site_coord[1] == 0: - base_position[0] = round(center_x, 1) - base_position[1] = round(center_y, 1) - else: - base_position[0] = round(site_coord[0], 1) - base_position[1] = round(site_coord[1], 1) - - # 打印检查 - print(f"min_z (adsorbate): {min_z}, max_z (slab): {max_z}, new base_position: {base_position}") - - # 平移吸附分子 + base_position[2] = max_z - min_z + 2 # 吸附物种最低原子处于slab以上2.2A，随后再尝试降低 a.translate(base_position) # 将分子平移到指定的 site_coord 位置 @@ -771,8 +814,7 @@ def _single_adsorption( # Add graph connections for metal_index in self.index[u]: slabs_list[-1].graph.add_edge(metal_index, bond + n) - ''' - # 评估当前构型并不断尝试将吸附物降低以尽可能贴近表面，直到score不再提高 + '''# 评估当前构型并不断尝试将吸附物降低以尽可能贴近表面，直到score不再提高 score_tmp = utils.score_configuration_hetero(coords=slabs_list[-1].get_positions(), symbols=slabs_list[-1].get_chemical_symbols(), z_surf=max_z) @@ -798,8 +840,8 @@ def _single_adsorption( str_log = str(ia) + ' | score = ' + str(np.round(score_configurations[-1],3)).ljust(8) + '(x,y,z): ' + str(base_position) ### print(str_log) with open('score_log.txt', 'a') as f: - f.write(str_log+"\n") - ''' + f.write(str_log+'\n') + ''' with open('score_log.txt', 'a') as f: f.write('\n') f.write('Ranking configurations by their scores:\n') @@ -809,8 +851,9 @@ def _single_adsorption( return slabs_list #lbx - '''if base_position[2] == 0.0: - z_coordinates = rotated_positions[:, 2] + if base_position[2] == 0.0: + #z_coordinates = rotated_positions[:, 2] + z_coordinates = atoms.get_positions()[:, 2] min_z = np.min(z_coordinates) #得到分子z轴最小值 #print(rotated_positions) #print(min_z) @@ -840,7 +883,7 @@ def _single_adsorption( # Add graph connections for metal_index in self.index[u]: - slab.graph.add_edge(metal_index, bond + n)''' + slab.graph.add_edge(metal_index, bond + n) return slab From f2bb1c2813d4e5a6ef302d6de521fe7ff024736a Mon Sep 17 00:00:00 2001 From: zqxian <2821838490@qq.com> Date: Sun, 22 Jun 2025 22:19:10 +0800 Subject: [PATCH 2/2] =?UTF-8?q?=E4=BF=9D=E7=95=99=E4=BD=8D=E7=82=B9?= =?UTF-8?q?=E5=92=8C=E6=A8=A1=E5=BC=8F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- HTMACat/catkit/gen/adsorption.py | 131 ++++++------------------------- 1 file changed, 26 insertions(+), 105 deletions(-) diff --git a/HTMACat/catkit/gen/adsorption.py b/HTMACat/catkit/gen/adsorption.py index 08fedaf..8122660 100644 --- a/HTMACat/catkit/gen/adsorption.py +++ b/HTMACat/catkit/gen/adsorption.py @@ -1,12 +1,6 @@ from . import defaults from . import utils from . import symmetry -from ase.build import molecule, fcc111 -from ase.io import write -from ase.thermochemistry import IdealGasThermo -from ase.vibrations import Vibrations - -from copy import deepcopy import HTMACat.catkit as catkit import matplotlib.pyplot as plt import itertools @@ -573,22 +567,21 @@ def get_principal_axes(inertia_tensor): return principal_axes def _single_adsorption( - self, - adsorbate, - bond, - slab=None, - site_index=0, - auto_construct=True, - enable_rotate_xoy=True, - rotation_mode='vnn', - rotation_args={}, - direction_mode='default', - direction_args={}, # 用于指定惯性矩模式 - site_coord=None, - z_bias=0, - symmetric=True, - verbose=False - ): + self, + adsorbate, + bond, + slab=None, + site_index=0, + auto_construct=True, + enable_rotate_xoy=True, + rotation_mode='vertical to vec_to_neigh_imgsite', + rotation_args={}, + direction_mode='default', # wzj 20230524 指定确定位向的方式 + direction_args={}, # wzj 20230524 为后续扩展预留的参数 + site_coord=None, + z_bias=0, + symmetric=True, + verbose=False): """Bond and adsorbate by a single atom.""" if slab is None: slab = self.slab.copy() @@ -610,82 +603,21 @@ def _single_adsorption( numbers = atoms.numbers[bond] R = radii[numbers] base_position = utils.trilaterate(top_sites[u], r + R, vector) - + # Zhaojie Wang 20230910(precise adsorption coord) if not site_coord is None: base_position = site_coord - + + branches = nx.bfs_successors(atoms.graph, bond) + atoms.translate(-atoms.positions[bond]) + + # Zhaojie Wang 20230510(direction), 20230828(rotation) + #lbx if auto_construct: - if direction_mode == 'asphericity': - # 根据分子的形状调整朝向，使其“平躺”在表面上 - masses = atoms.get_masses() - positions = atoms.get_positions() - #xzq - # 计算惯性矩张量 - inertia_tensor_value = self.inertia_tensor(positions, masses) - - # 获取主惯性轴 - principal_axes = self.get_principal_axes(inertia_tensor_value) - - # 自动处理三种惯性矩方向 - slabs_list = [] - for inertia_mode in range(1, 4): # 遍历第一、第二、第三惯性矩 - atoms_copy = atoms.copy() # ✅ 每次都从原始结构复制 - base_position = [0.0, 0.0, 0.0] - adsorption_vector = principal_axes[:, inertia_mode - 1] - atoms_copy.rotate(adsorption_vector, [0, 0, 1]) - - if enable_rotate_xoy and rotation_mode == 'vnn' and rotation_args != {}: - principle_axe = utils.solve_principle_axe_pca(atoms_copy.get_positions()) - if abs(rotation_args['vec_to_neigh_imgsite'][0]) < 1e-8: - target_vec = [1, 0, 0] - elif abs(rotation_args['vec_to_neigh_imgsite'][1]) < 1e-8: - target_vec = [0, 1, 0] - else: - target_vec = [-1/rotation_args['vec_to_neigh_imgsite'][0], - 1/rotation_args['vec_to_neigh_imgsite'][1], 0] - atoms_copy.rotate([principle_axe[0], principle_axe[1], 0], target_vec) - # =============== ✅ 关键部分：设置 Z 高度 =============== - z_coordinates = atoms_copy.get_positions()[:, 2] - min_z = np.min(z_coordinates) - - final_positions = slab.get_positions() - max_z = np.max(final_positions[:, 2]) - - - if abs(site_coord[2]) < 1e-6: # 如果为 0 或非常接近 0 - effective_distance = 2.5 - else: - effective_distance = site_coord[2] - - base_position[2] = round(max_z - min_z + effective_distance, 2) - - # xzq如果输入的 xy 为 [0, 0]，则将分子放在 slab 的 xy 中心 - if site_coord[0] == 0 and site_coord[1] == 0: - center_x, center_y = utils.center_slab(final_positions) # 计算 slab 的 xy 中心 - base_position[0] = round(center_x, 1) # 将分子放置在 xy 中心 - base_position[1] = round(center_y, 1) - else: - # 否则，使用输入的 site_coord 作为 xy 坐标 - base_position[0] = round(site_coord[0], 1) - base_position[1] = round(site_coord[1], 1) - - atoms_copy.translate(base_position) # 确保所有构型都以正确的 z 轴间距平移 - - # 将分子平移到指定的 site_coord 位置 - vec_tmp = site_coord - atoms_copy.get_positions()[bond] # 计算平移量 - atoms_copy.translate([vec_tmp[0], vec_tmp[1], 0]) # 只平移 x-y，不改变 z - - - slab_with_adsorbate = slab.copy() - slab_with_adsorbate += atoms_copy - slabs_list.append(slab_with_adsorbate) - - return slabs_list # 返回三种吸附构型的列表 - - elif direction_mode == 'bond_atom': + if direction_mode == 'bond_atom': # 根据参与吸附的原子确定位向，将物种“扶正” adsorption_vector, flag = utils.solve_normal_vector_linearsvc(atoms.get_positions(), bond) + ### print('adsorption_vector:\n', adsorption_vector) atoms.rotate(adsorption_vector, [0, 0, 1]) elif direction_mode == 'asphericity': # 根据分子的形状调整朝向，使其“平躺”在表面上 @@ -796,20 +728,12 @@ def _single_adsorption( n = len(slab) slabs_list = [] score_configurations = [] # 各个吸附构型(slab)的“分数” - for ia, a_orig in enumerate(atoms_list): + for ia,a in enumerate(atoms_list): slabs_list.append(copy.deepcopy(slab)) - a = copy.deepcopy(a_orig) # 确保每次处理原始未改动的 adsorbate dealt_positions = a.get_positions() min_z = np.min(dealt_positions[:,2]) #得到分子z轴最小值 base_position[2] = max_z - min_z + 2 # 吸附物种最低原子处于slab以上2.2A，随后再尝试降低 a.translate(base_position) - - # 将分子平移到指定的 site_coord 位置 - xy_translation = site_coord[:2] - a.get_positions()[bond][:2] - a.translate([xy_translation[0], xy_translation[1], 0]) # 只移动 x-y，不改变 z - print('Final base_position:', base_position) - - # 组合 slab 和吸附物 slabs_list[-1] += a # Add graph connections for metal_index in self.index[u]: @@ -847,7 +771,6 @@ def _single_adsorption( f.write('Ranking configurations by their scores:\n') for i,idx in enumerate(np.argsort(score_configurations)[::-1]): f.write(str(i).ljust(4)+': '+str(idx).ljust(8)+str(score_configurations[idx])+'\n') - return slabs_list #lbx @@ -860,8 +783,7 @@ def _single_adsorption( final_positions = slab.get_positions() #slab坐标 z_coordinates = final_positions[:, 2] max_z = np.max(z_coordinates) #获取slabz轴最大值 - base_position[2] = round(0 - min_z + 2.5 + max_z,1) - print(f"min_z (adsorbate): {min_z}, max_z (slab): {max_z}, new base_position[2]: {base_position[2]}") + base_position[2] = round(0 - min_z + 4.0 + max_z,1) #计算slab中心坐标 center_x, center_y = utils.center_slab(final_positions) #print("(x, y):", center_x,center_y,base_position[2]) @@ -887,7 +809,6 @@ def _single_adsorption( return slab - def _double_adsorption(self, adsorbate, bonds=None, edge_index=0): """Bond and adsorbate by two adjacent atoms.""" slab = self.slab.copy()