diff --git a/HTMACat/catkit/gen/adsorption.py b/HTMACat/catkit/gen/adsorption.py index 14792e9..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 @@ -571,23 +565,23 @@ 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, - 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() @@ -609,13 +603,23 @@ 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': + 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': # 根据分子的形状调整朝向，使其“平躺”在表面上 masses = atoms.get_masses() positions = atoms.get_positions() @@ -644,7 +648,7 @@ def _single_adsorption( 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) @@ -653,21 +657,15 @@ def _single_adsorption( if abs(site_coord[2]) < 1e-6: # 如果为 0 或非常接近 0 - effective_distance = 2.5 + effective_distance = 2 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) + #使用输入的 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 轴间距平移 @@ -680,13 +678,7 @@ def _single_adsorption( slab_with_adsorbate += atoms_copy slabs_list.append(slab_with_adsorbate) - return slabs_list # 返回三种吸附构型的列表 - - elif direction_mode == 'bond_atom': - # 根据参与吸附的原子确定位向，将物种“扶正” - adsorption_vector, flag = utils.solve_normal_vector_linearsvc(atoms.get_positions(), bond) - atoms.rotate(adsorption_vector, [0, 0, 1]) - + return slabs_list # 返回三种吸附构型的列表) elif direction_mode == 'hetero': # zjwang 20240815 # 适用于有明确“官能团”的偏链状的分子 # 求出从分子杂原子中心到所有原子中心的吸附矢量vec_p，将vec_p旋转到[0,0,1]竖直向上 @@ -736,43 +728,17 @@ 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轴最小值 - # 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 位置 - 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]: 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,27 +764,26 @@ 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') 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 - '''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) 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]) @@ -840,11 +805,10 @@ 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 - def _double_adsorption(self, adsorbate, bonds=None, edge_index=0): """Bond and adsorbate by two adjacent atoms.""" slab = self.slab.copy()