Skip to content
Merged

Dev #90

Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
126 changes: 45 additions & 81 deletions HTMACat/catkit/gen/adsorption.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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()
Expand All @@ -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()
Expand Down Expand Up @@ -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)

Expand All @@ -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 轴间距平移

Expand All @@ -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]竖直向上
Expand Down Expand Up @@ -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)
Expand All @@ -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])
Expand All @@ -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()
Expand Down