import numpy as np
from pymatgen.io.vasp import Vasprun
from tqdm import tqdm
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import os
import matplotlib as mpl
设置全局绘图参数
mpl.rcParams[‘font.family’] = ‘serif’ # 期刊常用字体
mpl.rcParams[‘mathtext.fontset’] = ‘dejavuserif’ # 数学字体
mpl.rcParams[‘axes.linewidth’] = 1.5 # 增加轴线宽度
mpl.rcParams[‘lines.linewidth’] = 2.5 # 线宽
mpl.rcParams[‘figure.dpi’] = 300 # 高分辨率
def get_angle(A, B, C):
“”“计算三个点之间的角度（B为顶点）”“”
BA = A - B
BC = C - B
cos_theta = np.dot(BA, BC) / (np.linalg.norm(BA) * np.linalg.norm(BC))
cos_theta = np.clip(cos_theta, -1.0, 1.0) # 避免数值误差
return np.degrees(np.arccos(cos_theta))
def analyze_system(system_name, vasprun_file, output_dir):
“”“分析单个体系的氢键分布”“”
# 创建体系特定输出目录
system_dir = os.path.join(output_dir, system_name)
os.makedirs(system_dir, exist_ok=True)
# 初始化分类存储 bond_types = { 'O-O': {'angles': [], 'color': '#1f77b4'}, 'O-F': {'angles': [], 'color': '#ff7f0e'}, 'F-O': {'angles': [], 'color': '#2ca02c'}, 'F-F': {'angles': [], 'color': '#d62728'} } # 参数设置 (可根据不同体系调整) step = 100 # 采样步长 donor_cutoff = 1.2 # 供体原子最大距离 (Å) receptor_cutoff = 2.5 # 受体原子最大距离 (Å) angle_threshold = 120 # 最小氢键角度 (度) print(f"[{system_name}] 正在读取 {vasprun_file} 文件...") try: vr = Vasprun(vasprun_file) except Exception as e: print(f"[{system_name}] 读取vasprun文件时出错: {e}") return system_name, bond_types # 返回空结果 structures = vr.structures print(f"[{system_name}] 总帧数: {len(structures)}, 采样步长: {step}") # 主分析循环 for frame_idx, structure in enumerate(tqdm(structures[::step], desc=f"{system_name}氢键分析", unit="帧")): # 预处理原子索引（优化性能） o_f_sites = [site for site in structure if site.species_string in ("O", "F")] for H_site in structure: if H_site.species_string != "H": continue # 寻找供体原子 donor = min( (site for site in o_f_sites), key=lambda x: H_site.distance(x), default=None ) if donor is None or H_site.distance(donor) > donor_cutoff: continue # 获取受体原子 receptors = [] for neighbor in structure.get_neighbors(H_site, r=receptor_cutoff): try: site = neighbor.site dist = neighbor.distance except AttributeError: site, dist, _, _ = neighbor if (site.species_string in ("O", "F") and site != donor and dist <= receptor_cutoff): receptors.append(site) # 分类计算角度 for receptor in receptors: angle = get_angle(donor.coords, H_site.coords, receptor.coords) if angle < angle_threshold: continue # 确定键类型 donor_type = donor.species_string receptor_type = receptor.species_string bond_key = f"{donor_type}-{receptor_type}" if bond_key in bond_types: bond_types[bond_key]['angles'].append(angle) # 保存体系统计数据 stats_data = [] for bond_type, data in bond_types.items(): angles = data['angles'] if len(angles) == 0: continue stats_data.append({ 'Type': bond_type, 'Count': len(angles), 'Mean': np.mean(angles), 'Std': np.std(angles), 'Median': np.median(angles), 'Min': np.min(angles), 'Max': np.max(angles) }) stats_df = pd.DataFrame(stats_data) stats_file = os.path.join(system_dir, "bond_type_stats.csv") stats_df.to_csv(stats_file, index=False) print(f"[{system_name}] 键型统计已保存到 {stats_file}") # 保存当前体系的分布图 plot_single_system(system_name, bond_types, system_dir) return system_name, bond_types
def plot_single_system(system_name, bond_types, output_dir):
“”“为单个体系绘制分布图并保存”“”
plt.figure(figsize=(10, 6))
sns.set_style(“whitegrid”)
# 统一设置颜色映射 bond_colors = { 'O-O': '#1f77b4', 'O-F': '#ff7f0e', 'F-O': '#2ca02c', 'F-F': '#d62728' } bins = np.linspace(90, 180, 30) # 绘制分布 for bond_type, data in bond_types.items(): angles = data['angles'] if len(angles) == 0: continue # KDE曲线 sns.kdeplot(angles, label=f"{bond_type} (n={len(angles)})", color=bond_colors[bond_type], linewidth=2.5, alpha=0.8) # 直方图叠加 plt.hist(angles, bins=bins, color=bond_colors[bond_type], alpha=0.2, density=True) # 设置图表属性 plt.xlabel("Hydrogen Bond Angle (degrees)", fontsize=14) plt.ylabel("Probability Density", fontsize=14) plt.title(f"Hydrogen Bond Angle Distribution: {system_name}", fontsize=16) plt.xlim(90, 180) plt.ylim(0, 0.035) plt.legend(title="Bond Type", fontsize=12) # 保存图像 plot_file = os.path.join(output_dir, f"{system_name}_bond_distribution.png") plt.tight_layout() plt.savefig(plot_file, dpi=300, bbox_inches='tight') plt.close() # 关闭图形释放内存 print(f"[{system_name}] 已保存分布图到 {plot_file}")
def plot_combined_results(results, output_dir):
“”“将四个体系的结果绘制在统一坐标的子图中并保存”“”
plt.figure(figsize=(14, 12))
sns.set(style=“whitegrid”, font_scale=1.2, palette=“muted”)
bins = np.linspace(90, 180, 30)
# 创建2x2子图 fig, axes = plt.subplots(2, 2, figsize=(14, 12), sharex=True, sharey=True) axes = axes.flatten() # 统一设置颜色映射 bond_colors = { 'O-O': '#1f77b4', 'O-F': '#ff7f0e', 'F-O': '#2ca02c', 'F-F': '#d62728' } # 为每个体系绘制子图 for idx, (system_name, bond_types) in enumerate(results.items()): ax = axes[idx] # 绘制每种键型的分布 for bond_type, data in bond_types.items(): angles = data['angles'] if len(angles) == 0: continue # KDE曲线 sns.kdeplot(angles, label=f"{bond_type}", color=bond_colors[bond_type], linewidth=2.5, ax=ax) # 直方图叠加 ax.hist(angles, bins=bins, color=bond_colors[bond_type], alpha=0.2, density=True) # 设置子图标题和标签 ax.set_title(f"{system_name}", fontsize=14, fontweight='bold') ax.set_xlabel("Hydrogen Bond Angle (degrees)", fontsize=12) ax.set_ylabel("Probability Density", fontsize=12) ax.set_xlim(90, 180) ax.set_ylim(0, 0.035) # 统一Y轴范围 ax.legend(title="Bond Type", fontsize=10, loc='upper left') ax.grid(True, linestyle='--', alpha=0.7) # 添加统计信息文本框 stats_text = "\n".join([f"{bt}: n={len(data['angles'])}" for bt, data in bond_types.items() if len(data['angles']) > 0]) ax.text(0.95, 0.95, stats_text, transform=ax.transAxes, verticalalignment='top', horizontalalignment='right', bbox=dict(boxstyle='round', facecolor='white', alpha=0.8), fontsize=9) # 添加整体标题 plt.suptitle("Hydrogen Bond Angle Distribution in Different Systems", fontsize=16, fontweight='bold', y=0.98) # 调整布局 plt.tight_layout(rect=[0, 0, 1, 0.96]) # 保存组合图 combined_file = os.path.join(output_dir, "combined_bond_distribution.png") plt.savefig(combined_file, dpi=300, bbox_inches='tight') plt.close() # 关闭图形释放内存 print(f"已保存组合分布图到 {combined_file}") # 保存组合图PDF版本（期刊投稿要求） combined_pdf = os.path.join(output_dir, "combined_bond_distribution.pdf") plt.savefig(combined_pdf, bbox_inches='tight') print(f"已保存组合分布图PDF版本到 {combined_pdf}")
if name == “main”:
# ===== 配置四个体系 =====
# 假设所有vasprun文件都在当前目录下
systems = {
“System_1”: “vasprun1.xml”,
“System_2”: “vasprun2.xml”,
“System_3”: “vasprun3.xml”,
“System_4”: “vasprun4.xml”
}
output_dir = "hydrogen_bond_analysis" os.makedirs(output_dir, exist_ok=True) # ===== 顺序处理四个体系 ===== results = {} for system_name, file in tqdm(systems.items(), desc="处理体系"): try: system_name, bond_types = analyze_system(system_name, file, output_dir) results[system_name] = bond_types except Exception as e: print(f"处理体系 {system_name} 时出错: {e}") # ===== 绘制组合图 ===== plot_combined_results(results, output_dir) print("\n所有体系分析完成！结果保存在", output_dir) print("包含:") print(" - 每个体系的单独目录（包含统计数据和分布图）") print(" - 组合分布图 (combined_bond_distribution.png 和 .pdf)")以上代码实现了对四个vasprun体系中计算四类氢键键角的角度分布情况，在这里我们对其中的原子分类逻辑做进一步优化，对计算的类别进一步优化，然后将该代码转为在58核100g内存的台式ananconda promt中执行，先识别P，在P周围搜寻O原子，如果该O原子距离在1.6埃以内则视为Op，而对于Op在其周围搜寻H原子，如果该H在距离1.3埃以内则视为成键即该H为Hp，通过是否有Hp则可以识别P-OH与P=O/P-O，在这里P=O和P-O还不能区分，在这里我们将不对P=O进行区分识别，为方便清晰，将两者统一记为P-O/P=O。接着体系中全部的O原子在去除Op之后剩下的O，在这些剩余的O周围搜寻整体的H，如果H的距离在1.2埃以内则视为成键，然后依照成键的H数量判定：如果H的数量为1，则记为-OH羟基（在这里不需要计算羟基部分，只是识别出来有利于逻辑完整性，并不参与RDF计算，也不需要特别标注表明），H的数量为2，则记为H2O水（该O也随之记为Ow，对应的两个H也记为Hw），如果H的数量为3，则记为水合氢离子（该O随之记为Oh，对应的三个H也记为Hh）。体系中存在质子转移的情况，所以需要每一帧重新识别原子的归属问题，如果H同时处于两个成键识别范围则按照就近原则，离哪个近则归属到哪一个（这里包括磷酸-磷酸，磷酸-水，磷酸-水合氢离子，水-水，水-水合氢离子，水合氢离子-水合氢离子，如果H同时处于某种情况下两个化学成键范围则采用就近原则），在实时重新归属质子的情况下，计算出包含质子转移部分的RDF，在这里，我们将排除自身化学键的阈值先设置为0，不需要只看氢键部分了。直接将-OH视为质子转移或者不完整而直接忽略即可，磷酸上的O需要通过H来进一步识别，所以符合距离的氧可暂时作为Op的候选，等H的识别完成再进行细分P=O/P-O和P-OH。同时由于后续RDF的计算过程中我们发现，在这里我们只需要将原先对四种类别的计算改为对九类的计算，分别为P-O/P=O···Hw P-O/P=O···Hh P-O/P=O···Hp P-OH···Ow Hw···Ow Hh···Ow Hw···F Hh···F Hp···F，以上补充的逻辑中提供了化学键的判别（供体和H之间的距离），在这里我们对H和受体之间的距离再进一步补充，这九类分别的距离分别为2.375 2.275 2.175 2.275 2.375 2.225 2.225 2.175 2.275.这数据的导出需要能导入origin中方便我后续作图，在原代码中文件名后加上2，作为第二版的意思。
回答
向我提问的人太多了。正在努力扩容中,请稍后再试。
import numpy as np
from pymatgen.io.vasp import Vasprun
from tqdm import tqdm
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import os
import matplotlib as mpl
设置全局绘图参数
mpl.rcParams[‘font.family’] = ‘serif’ # 期刊常用字体
mpl.rcParams[‘mathtext.fontset’] = ‘dejavuserif’ # 数学字体
mpl.rcParams[‘axes.linewidth’] = 1.5 # 增加轴线宽度
mpl.rcParams[‘lines.linewidth’] = 2.5 # 线宽
mpl.rcParams[‘figure.dpi’] = 300 # 高分辨率
def get_angle(A, B, C):
“”“计算三个点之间的角度（B为顶点）”“”
BA = A - B
BC = C - B
cos_theta = np.dot(BA, BC) / (np.linalg.norm(BA) * np.linalg.norm(BC))
cos_theta = np.clip(cos_theta, -1.0, 1.0) # 避免数值误差
return np.degrees(np.arccos(cos_theta))
def analyze_system(system_name, vasprun_file, output_dir):
“”“分析单个体系的氢键分布”“”
创建体系特定输出目录
system_dir = os.path.join(output_dir, system_name)
os.makedirs(system_dir, exist_ok=True)
初始化分类存储 bond_types = { ‘O-O’: {‘angles’: [], ‘color’: ‘#1f77b4’}, ‘O-F’: {‘angles’: [], ‘color’: ‘#ff7f0e’}, ‘F-O’: {‘angles’: [], ‘color’: ‘#2ca02c’}, ‘F-F’: {‘angles’: [], ‘color’: ‘#d62728’} } # 参数设置 (可根据不同体系调整) step = 100 # 采样步长 donor_cutoff = 1.2 # 供体原子最大距离 (Å) receptor_cutoff = 2.5 # 受体原子最大距离 (Å) angle_threshold = 120 # 最小氢键角度 (度) print(f"[{system_name}] 正在读取 {vasprun_file} 文件…“) try: vr = Vasprun(vasprun_file) except Exception as e: print(f”[{system_name}] 读取vasprun文件时出错: {e}“) return system_name, bond_types # 返回空结果 structures = vr.structures print(f”[{system_name}] 总帧数: {len(structures)}, 采样步长: {step}“) # 主分析循环 for frame_idx, structure in enumerate(tqdm(structures[::step], desc=f”{system_name}氢键分析", unit=“帧”)): # 预处理原子索引（优化性能） o_f_sites = [site for site in structure if site.species_string in (“O”, “F”)] for H_site in structure: if H_site.species_string != “H”: continue # 寻找供体原子 donor = min( (site for site in o_f_sites), key=lambda x: H_site.distance(x), default=None ) if donor is None or H_site.distance(donor) > donor_cutoff: continue # 获取受体原子 receptors = [] for neighbor in structure.get_neighbors(H_site, r=receptor_cutoff): try: site = neighbor.site dist = neighbor.distance except AttributeError: site, dist, _, _ = neighbor if (site.species_string in (“O”, “F”) and site != donor and dist <= receptor_cutoff): receptors.append(site) # 分类计算角度 for receptor in receptors: angle = get_angle(donor.coords, H_site.coords, receptor.coords) if angle < angle_threshold: continue # 确定键类型 donor_type = donor.species_string receptor_type = receptor.species_string bond_key = f"{donor_type}-{receptor_type}" if bond_key in bond_types: bond_types[bond_key][‘angles’].append(angle) # 保存体系统计数据 stats_data = [] for bond_type, data in bond_types.items(): angles = data[‘angles’] if len(angles) == 0: continue stats_data.append({ ‘Type’: bond_type, ‘Count’: len(angles), ‘Mean’: np.mean(angles), ‘Std’: np.std(angles), ‘Median’: np.median(angles), ‘Min’: np.min(angles), ‘Max’: np.max(angles) }) stats_df = pd.DataFrame(stats_data) stats_file = os.path.join(system_dir, “bond_type_stats.csv”) stats_df.to_csv(stats_file, index=False) print(f"[{system_name}] 键型统计已保存到 {stats_file}") # 保存当前体系的分布图 plot_single_system(system_name, bond_types, system_dir) return system_name, bond_types
def plot_single_system(system_name, bond_types, output_dir):
“”“为单个体系绘制分布图并保存”“”
plt.figure(figsize=(10, 6))
sns.set_style(“whitegrid”)
统一设置颜色映射 bond_colors = { ‘O-O’: ‘#1f77b4’, ‘O-F’: ‘#ff7f0e’, ‘F-O’: ‘#2ca02c’, ‘F-F’: ‘#d62728’ } bins = np.linspace(90, 180, 30) # 绘制分布 for bond_type, data in bond_types.items(): angles = data[‘angles’] if len(angles) == 0: continue # KDE曲线 sns.kdeplot(angles, label=f"{bond_type} (n={len(angles)})“, color=bond_colors[bond_type], linewidth=2.5, alpha=0.8) # 直方图叠加 plt.hist(angles, bins=bins, color=bond_colors[bond_type], alpha=0.2, density=True) # 设置图表属性 plt.xlabel(“Hydrogen Bond Angle (degrees)”, fontsize=14) plt.ylabel(“Probability Density”, fontsize=14) plt.title(f"Hydrogen Bond Angle Distribution: {system_name}”, fontsize=16) plt.xlim(90, 180) plt.ylim(0, 0.035) plt.legend(title=“Bond Type”, fontsize=12) # 保存图像 plot_file = os.path.join(output_dir, f"{system_name}_bond_distribution.png") plt.tight_layout() plt.savefig(plot_file, dpi=300, bbox_inches=‘tight’) plt.close() # 关闭图形释放内存 print(f"[{system_name}] 已保存分布图到 {plot_file}")
def plot_combined_results(results, output_dir):
“”“将四个体系的结果绘制在统一坐标的子图中并保存”“”
plt.figure(figsize=(14, 12))
sns.set(style=“whitegrid”, font_scale=1.2, palette=“muted”)
bins = np.linspace(90, 180, 30)
创建2x2子图 fig, axes = plt.subplots(2, 2, figsize=(14, 12), sharex=True, sharey=True) axes = axes.flatten() # 统一设置颜色映射 bond_colors = { ‘O-O’: ‘#1f77b4’, ‘O-F’: ‘#ff7f0e’, ‘F-O’: ‘#2ca02c’, ‘F-F’: ‘#d62728’ } # 为每个体系绘制子图 for idx, (system_name, bond_types) in enumerate(results.items()): ax = axes[idx] # 绘制每种键型的分布 for bond_type, data in bond_types.items(): angles = data[‘angles’] if len(angles) == 0: continue # KDE曲线 sns.kdeplot(angles, label=f"{bond_type}“, color=bond_colors[bond_type], linewidth=2.5, ax=ax) # 直方图叠加 ax.hist(angles, bins=bins, color=bond_colors[bond_type], alpha=0.2, density=True) # 设置子图标题和标签 ax.set_title(f”{system_name}“, fontsize=14, fontweight=‘bold’) ax.set_xlabel(“Hydrogen Bond Angle (degrees)”, fontsize=12) ax.set_ylabel(“Probability Density”, fontsize=12) ax.set_xlim(90, 180) ax.set_ylim(0, 0.035) # 统一Y轴范围 ax.legend(title=“Bond Type”, fontsize=10, loc=‘upper left’) ax.grid(True, linestyle=‘–’, alpha=0.7) # 添加统计信息文本框 stats_text = “\n”.join([f”{bt}: n={len(data[‘angles’])}" for bt, data in bond_types.items() if len(data[‘angles’]) > 0]) ax.text(0.95, 0.95, stats_text, transform=ax.transAxes, verticalalignment=‘top’, horizontalalignment=‘right’, bbox=dict(boxstyle=‘round’, facecolor=‘white’, alpha=0.8), fontsize=9) # 添加整体标题 plt.suptitle(“Hydrogen Bond Angle Distribution in Different Systems”, fontsize=16, fontweight=‘bold’, y=0.98) # 调整布局 plt.tight_layout(rect=[0, 0, 1, 0.96]) # 保存组合图 combined_file = os.path.join(output_dir, “combined_bond_distribution.png”) plt.savefig(combined_file, dpi=300, bbox_inches=‘tight’) plt.close() # 关闭图形释放内存 print(f"已保存组合分布图到 {combined_file}“) # 保存组合图PDF版本（期刊投稿要求） combined_pdf = os.path.join(output_dir, “combined_bond_distribution.pdf”) plt.savefig(combined_pdf, bbox_inches=‘tight’) print(f"已保存组合分布图PDF版本到 {combined_pdf}”)
if name == “main”:
===== 配置四个体系 =====
假设所有vasprun文件都在当前目录下
systems = {
“System_1”: “vasprun1.xml”,
“System_2”: “vasprun2.xml”,
“System_3”: “vasprun3.xml”,
“System_4”: “vasprun4.xml”
}
output_dir = “hydrogen_bond_analysis” os.makedirs(output_dir, exist_ok=True) # ===== 顺序处理四个体系 ===== results = {} for system_name, file in tqdm(systems.items(), desc=“处理体系”): try: system_name, bond_types = analyze_system(system_name, file, output_dir) results[system_name] = bond_types except Exception as e: print(f"处理体系 {system_name} 时出错: {e}“) # ===== 绘制组合图 ===== plot_combined_results(results, output_dir) print(”\n所有体系分析完成！结果保存在", output_dir) print(“包含:”) print(" - 每个体系的单独目录（包含统计数据和分布图）“) print(” - 组合分布图 (combined_bond_distribution.png 和 .pdf)")以上代码实现了对四个vasprun体系中计算四类氢键键角的角度分布情况，在这里我们对其中的原子分类逻辑做进一步优化，对计算的类别进一步优化，然后将该代码转为在58核100g内存的台式ananconda promt中执行，先识别P，在P周围搜寻O原子，如果该O原子距离在1.6埃以内则视为Op，而对于Op在其周围搜寻H原子，如果该H在距离1.3埃以内则视为成键即该H为Hp，通过是否有Hp则可以识别P-OH与P=O/P-O，在这里P=O和P-O还不能区分，在这里我们将不对P=O进行区分识别，为方便清晰，将两者统一记为P-O/P=O。接着体系中全部的O原子在去除Op之后剩下的O，在这些剩余的O周围搜寻整体的H，如果H的距离在1.2埃以内则视为成键，然后依照成键的H数量判定：如果H的数量为1，则记为-OH羟基（在这里不需要计算羟基部分，只是识别出来有利于逻辑完整性，并不参与RDF计算，也不需要特别标注表明），H的数量为2，则记为H2O水（该O也随之记为Ow，对应的两个H也记为Hw），如果H的数量为3，则记为水合氢离子（该O随之记为Oh，对应的三个H也记为Hh）。体系中存在质子转移的情况，所以需要每一帧重新识别原子的归属问题，如果H同时处于两个成键识别范围则按照就近原则，离哪个近则归属到哪一个（这里包括磷酸-磷酸，磷酸-水，磷酸-水合氢离子，水-水，水-水合氢离子，水合氢离子-水合氢离子，如果H同时处于某种情况下两个化学成键范围则采用就近原则），在实时重新归属质子的情况下，计算出包含质子转移部分的RDF，在这里，我们将排除自身化学键的阈值先设置为0，不需要只看氢键部分了。直接将-OH视为质子转移或者不完整而直接忽略即可，磷酸上的O需要通过H来进一步识别，所以符合距离的氧可暂时作为Op的候选，等H的识别完成再进行细分P=O/P-O和P-OH。同时由于后续RDF的计算过程中我们发现，在这里我们只需要将原先对四种类别的计算改为对九类的计算，分别为P-O/P=O···Hw P-O/P=O···Hh P-O/P=O···Hp P-OH···Ow Hw···Ow Hh···Ow Hw···F Hh···F Hp···F，以上补充的逻辑中提供了化学键的判别（供体和H之间的距离），在这里我们对H和受体之间的距离再进一步补充，这九类分别的距离分别为2.375 2.275 2.175 2.275 2.375 2.225 2.225 2.175 2.275.这数据的导出需要能导入origin中方便我后续作图，在原代码中文件名后加上2，作为第二版的意思。
回答
向我提问的人太多了。正在努力扩容中,请稍后再试。
import numpy as np
from pymatgen.io.vasp import Vasprun
from tqdm import tqdm
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import os
import matplotlib as mpl
设置全局绘图参数
mpl.rcParams[‘font.family’] = ‘serif’ # 期刊常用字体
mpl.rcParams[‘mathtext.fontset’] = ‘dejavuserif’ # 数学字体
mpl.rcParams[‘axes.linewidth’] = 1.5 # 增加轴线宽度
mpl.rcParams[‘lines.linewidth’] = 2.5 # 线宽
mpl.rcParams[‘figure.dpi’] = 300 # 高分辨率
def get_angle(A, B, C):
“”“计算三个点之间的角度（B为顶点）”“”
BA = A - B
BC = C - B
cos_theta = np.dot(BA, BC) / (np.linalg.norm(BA) * np.linalg.norm(BC))
cos_theta = np.clip(cos_theta, -1.0, 1.0) # 避免数值误差
return np.degrees(np.arccos(cos_theta))
def analyze_system(system_name, vasprun_file, output_dir):
“”“分析单个体系的氢键分布”“”
创建体系特定输出目录
system_dir = os.path.join(output_dir, system_name)
os.makedirs(system_dir, exist_ok=True)
初始化分类存储 bond_types = { ‘O-O’: {‘angles’: [], ‘color’: ‘#1f77b4’}, ‘O-F’: {‘angles’: [], ‘color’: ‘#ff7f0e’}, ‘F-O’: {‘angles’: [], ‘color’: ‘#2ca02c’}, ‘F-F’: {‘angles’: [], ‘color’: ‘#d62728’} } # 参数设置 (可根据不同体系调整) step = 100 # 采样步长 donor_cutoff = 1.2 # 供体原子最大距离 (Å) receptor_cutoff = 2.5 # 受体原子最大距离 (Å) angle_threshold = 120 # 最小氢键角度 (度) print(f"[{system_name}] 正在读取 {vasprun_file} 文件…“) try: vr = Vasprun(vasprun_file) except Exception as e: print(f”[{system_name}] 读取vasprun文件时出错: {e}“) return system_name, bond_types # 返回空结果 structures = vr.structures print(f”[{system_name}] 总帧数: {len(structures)}, 采样步长: {step}“) # 主分析循环 for frame_idx, structure in enumerate(tqdm(structures[::step], desc=f”{system_name}氢键分析", unit=“帧”)): # 预处理原子索引（优化性能） o_f_sites = [site for site in structure if site.species_string in (“O”, “F”)] for H_site in structure: if H_site.species_string != “H”: continue # 寻找供体原子 donor = min( (site for site in o_f_sites), key=lambda x: H_site.distance(x), default=None ) if donor is None or H_site.distance(donor) > donor_cutoff: continue # 获取受体原子 receptors = [] for neighbor in structure.get_neighbors(H_site, r=receptor_cutoff): try: site = neighbor.site dist = neighbor.distance except AttributeError: site, dist, _, _ = neighbor if (site.species_string in (“O”, “F”) and site != donor and dist <= receptor_cutoff): receptors.append(site) # 分类计算角度 for receptor in receptors: angle = get_angle(donor.coords, H_site.coords, receptor.coords) if angle < angle_threshold: continue # 确定键类型 donor_type = donor.species_string receptor_type = receptor.species_string bond_key = f"{donor_type}-{receptor_type}" if bond_key in bond_types: bond_types[bond_key][‘angles’].append(angle) # 保存体系统计数据 stats_data = [] for bond_type, data in bond_types.items(): angles = data[‘angles’] if len(angles) == 0: continue stats_data.append({ ‘Type’: bond_type, ‘Count’: len(angles), ‘Mean’: np.mean(angles), ‘Std’: np.std(angles), ‘Median’: np.median(angles), ‘Min’: np.min(angles), ‘Max’: np.max(angles) }) stats_df = pd.DataFrame(stats_data) stats_file = os.path.join(system_dir, “bond_type_stats.csv”) stats_df.to_csv(stats_file, index=False) print(f"[{system_name}] 键型统计已保存到 {stats_file}") # 保存当前体系的分布图 plot_single_system(system_name, bond_types, system_dir) return system_name, bond_types
def plot_single_system(system_name, bond_types, output_dir):
“”“为单个体系绘制分布图并保存”“”
plt.figure(figsize=(10, 6))
sns.set_style(“whitegrid”)
统一设置颜色映射 bond_colors = { ‘O-O’: ‘#1f77b4’, ‘O-F’: ‘#ff7f0e’, ‘F-O’: ‘#2ca02c’, ‘F-F’: ‘#d62728’ } bins = np.linspace(90, 180, 30) # 绘制分布 for bond_type, data in bond_types.items(): angles = data[‘angles’] if len(angles) == 0: continue # KDE曲线 sns.kdeplot(angles, label=f"{bond_type} (n={len(angles)})“, color=bond_colors[bond_type], linewidth=2.5, alpha=0.8) # 直方图叠加 plt.hist(angles, bins=bins, color=bond_colors[bond_type], alpha=0.2, density=True) # 设置图表属性 plt.xlabel(“Hydrogen Bond Angle (degrees)”, fontsize=14) plt.ylabel(“Probability Density”, fontsize=14) plt.title(f"Hydrogen Bond Angle Distribution: {system_name}”, fontsize=16) plt.xlim(90, 180) plt.ylim(0, 0.035) plt.legend(title=“Bond Type”, fontsize=12) # 保存图像 plot_file = os.path.join(output_dir, f"{system_name}_bond_distribution.png") plt.tight_layout() plt.savefig(plot_file, dpi=300, bbox_inches=‘tight’) plt.close() # 关闭图形释放内存 print(f"[{system_name}] 已保存分布图到 {plot_file}")
def plot_combined_results(results, output_dir):
“”“将四个体系的结果绘制在统一坐标的子图中并保存”“”
plt.figure(figsize=(14, 12))
sns.set(style=“whitegrid”, font_scale=1.2, palette=“muted”)
bins = np.linspace(90, 180, 30)
创建2x2子图 fig, axes = plt.subplots(2, 2, figsize=(14, 12), sharex=True, sharey=True) axes = axes.flatten() # 统一设置颜色映射 bond_colors = { ‘O-O’: ‘#1f77b4’, ‘O-F’: ‘#ff7f0e’, ‘F-O’: ‘#2ca02c’, ‘F-F’: ‘#d62728’ } # 为每个体系绘制子图 for idx, (system_name, bond_types) in enumerate(results.items()): ax = axes[idx] # 绘制每种键型的分布 for bond_type, data in bond_types.items(): angles = data[‘angles’] if len(angles) == 0: continue # KDE曲线 sns.kdeplot(angles, label=f"{bond_type}“, color=bond_colors[bond_type], linewidth=2.5, ax=ax) # 直方图叠加 ax.hist(angles, bins=bins, color=bond_colors[bond_type], alpha=0.2, density=True) # 设置子图标题和标签 ax.set_title(f”{system_name}“, fontsize=14, fontweight=‘bold’) ax.set_xlabel(“Hydrogen Bond Angle (degrees)”, fontsize=12) ax.set_ylabel(“Probability Density”, fontsize=12) ax.set_xlim(90, 180) ax.set_ylim(0, 0.035) # 统一Y轴范围 ax.legend(title=“Bond Type”, fontsize=10, loc=‘upper left’) ax.grid(True, linestyle=‘–’, alpha=0.7) # 添加统计信息文本框 stats_text = “\n”.join([f”{bt}: n={len(data[‘angles’])}" for bt, data in bond_types.items() if len(data[‘angles’]) > 0]) ax.text(0.95, 0.95, stats_text, transform=ax.transAxes, verticalalignment=‘top’, horizontalalignment=‘right’, bbox=dict(boxstyle=‘round’, facecolor=‘white’, alpha=0.8), fontsize=9) # 添加整体标题 plt.suptitle(“Hydrogen Bond Angle Distribution in Different Systems”, fontsize=16, fontweight=‘bold’, y=0.98) # 调整布局 plt.tight_layout(rect=[0, 0, 1, 0.96]) # 保存组合图 combined_file = os.path.join(output_dir, “combined_bond_distribution.png”) plt.savefig(combined_file, dpi=300, bbox_inches=‘tight’) plt.close() # 关闭图形释放内存 print(f"已保存组合分布图到 {combined_file}“) # 保存组合图PDF版本（期刊投稿要求） combined_pdf = os.path.join(output_dir, “combined_bond_distribution.pdf”) plt.savefig(combined_pdf, bbox_inches=‘tight’) print(f"已保存组合分布图PDF版本到 {combined_pdf}”)
if name == “main”:
===== 配置四个体系 =====
假设所有vasprun文件都在当前目录下
systems = {
“System_1”: “vasprun1.xml”,
“System_2”: “vasprun2.xml”,
“System_3”: “vasprun3.xml”,
“System_4”: “vasprun4.xml”
}
output_dir = “hydrogen_bond_analysis” os.makedirs(output_dir, exist_ok=True) # ===== 顺序处理四个体系 ===== results = {} for system_name, file in tqdm(systems.items(), desc=“处理体系”): try: system_name, bond_types = analyze_system(system_name, file, output_dir) results[system_name] = bond_types except Exception as e: print(f"处理体系 {system_name} 时出错: {e}“) # ===== 绘制组合图 ===== plot_combined_results(results, output_dir) print(”\n所有体系分析完成！结果保存在", output_dir) print(“包含:”) print(" - 每个体系的单独目录（包含统计数据和分布图）“) print(” - 组合分布图 (combined_bond_distribution.png 和 .pdf)")以上代码实现了对四个vasprun体系中计算四类氢键键角的角度分布情况，在这里我们对其中的原子分类逻辑做进一步优化，对计算的类别进一步优化，然后将该代码转为在58核100g内存的台式ananconda promt中执行，先识别P，在P周围搜寻O原子，如果该O原子距离在1.6埃以内则视为Op，而对于Op在其周围搜寻H原子，如果该H在距离1.3埃以内则视为成键即该H为Hp，通过是否有Hp则可以识别P-OH与P=O/P-O，在这里P=O和P-O还不能区分，在这里我们将不对P=O进行区分识别，为方便清晰，将两者统一记为P-O/P=O。接着体系中全部的O原子在去除Op之后剩下的O，在这些剩余的O周围搜寻整体的H，如果H的距离在1.2埃以内则视为成键，然后依照成键的H数量判定：如果H的数量为1，则记为-OH羟基（在这里不需要计算羟基部分，只是识别出来有利于逻辑完整性，并不参与RDF计算，也不需要特别标注表明），H的数量为2，则记为H2O水（该O也随之记为Ow，对应的两个H也记为Hw），如果H的数量为3，则记为水合氢离子（该O随之记为Oh，对应的三个H也记为Hh）。体系中存在质子转移的情况，所以需要每一帧重新识别原子的归属问题，如果H同时处于两个成键识别范围则按照就近原则，离哪个近则归属到哪一个（这里包括磷酸-磷酸，磷酸-水，磷酸-水合氢离子，水-水，水-水合氢离子，水合氢离子-水合氢离子，如果H同时处于某种情况下两个化学成键范围则采用就近原则），在实时重新归属质子的情况下，计算出包含质子转移部分的RDF，在这里，我们将排除自身化学键的阈值先设置为0，不需要只看氢键部分了。直接将-OH视为质子转移或者不完整而直接忽略即可，磷酸上的O需要通过H来进一步识别，所以符合距离的氧可暂时作为Op的候选，等H的识别完成再进行细分P=O/P-O和P-OH。同时由于后续RDF的计算过程中我们发现，在这里我们只需要将原先对四种类别的计算改为对九类的计算，分别为P-O/P=O···Hw P-O/P=O···Hh P-O/P=O···Hp P-OH···Ow Hw···Ow Hh···Ow Hw···F Hh···F Hp···F，以上补充的逻辑中提供了化学键的判别（供体和H之间的距离），在这里我们对H和受体之间的距离再进一步补充，这九类分别的距离分别为2.375 2.275 2.175 2.275 2.375 2.225 2.225 2.175 2.275.这数据的导出需要能导入origin中方便我后续作图，在原代码中文件名后加上2，作为第二版的意思。其中控制核数在58左右，内存控制在100g左右。
最新发布
本文介绍如何在R语言中自定义pairs函数以实现特定的数据可视化需求，包括修改面板显示、添加特殊标记等功能。
扫一扫
4566
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
