O_APPEND:原子操作

当多个进程打开同一个文件写入日志的时候,OPEN时指定了O_APPEND参数,UNIX能保证这个操作是原子的,程序不需要自己加锁

/*log1.c*/
#include<stdio.h>
#include<string.h>
#include<fcntl.h>
#include<unistd.h>

int Max = 300;
int main()
{
    int i;
    int fd;
    fd = open("1.log",O_WRONLY | O_APPEND);
    for(i=0;i<Max;i++)
    {
        char msg[100];
        sprintf(msg,"log1: message %d...\n",i);
        sleep(1);//让CPU调度走
        write(fd,msg,strlen(msg));
        write(STDOUT_FILENO,msg,strlen(msg));
    }
    close(fd);
    return 0;    
}

 

/*log2.c*/
#include<stdio.h>
#include<string.h>
#include<fcntl.h>
#include<unistd.h>

int Max = 300;
int main()
{
    int i;
    int fd;
    fd = open("1.log",O_WRONLY | O_APPEND);
    for(i=0;i<Max;i++)
    {
        char msg[100];
        sprintf(msg,"log2: message %d...\n",i);
        sleep(1);//让CPU调度走
        write(fd,msg,strlen(msg));
        write(STDOUT_FILENO,msg,strlen(msg));
    }
    close(fd);
    return 0;    
}

 

/*log3.c*/
#include<stdio.h>
#include<string.h>
#include<fcntl.h>
#include<unistd.h>

int Max = 300;
int main()
{
    int i;
    int fd;
    fd = open("2.log",O_WRONLY);
    for(i=0;i<Max;i++)
    {
        char msg[100];
        sprintf(msg,"log3: message %d...\n",i);
        lseek(fd,0,SEEK_END);
        sleep(1);//让CPU调度走
        write(fd,msg,strlen(msg));
        write(STDOUT_FILENO,msg,strlen(msg));
    }
    close(fd);
    return 0;    
}

 

/*log4.c*/
#include<stdio.h>
#include<string.h>
#include<fcntl.h>
#include<unistd.h>

int Max = 300;
int main()
{
    int i;
    int fd;
    fd = open("2.log",O_WRONLY);
    for(i=0;i<Max;i++)
    {
        char msg[100];
        sprintf(msg,"log4: message %d...\n",i);
        lseek(fd,0,SEEK_END);
        sleep(3);//让CPU调度走
        write(fd,msg,strlen(msg));
        write(STDOUT_FILENO,msg,strlen(msg));
    }
    close(fd);
    return 0;    
}

 

转载于:https://www.cnblogs.com/code-style/archive/2012/06/28/2566566.html

import numpy as np import matplotlib.pyplot as plt from pymatgen.io.vasp import Vasprun from pymatgen.core.structure import Structure from scipy.spatial import cKDTree from tqdm import tqdm import matplotlib as mpl import warnings import os import csv import argparse import multiprocessing from scipy.stats import gaussian_kde import time import sys # 忽略可能的警告 warnings.filterwarnings("ignore", category=UserWarning) # 专业绘图设置 plt.style.use('seaborn-v0_8-whitegrid') mpl.rcParams.update({ 'font.family': 'serif', 'font.serif': ['Times New Roman', 'DejaVu Serif'], 'font.size': 12, 'axes.labelsize': 14, 'axes.titlesize': 16, 'xtick.labelsize': 12, 'ytick.labelsize': 12, 'figure.dpi': 600, 'savefig.dpi': 600, 'figure.figsize': (8, 6), 'lines.linewidth': 2.0, 'legend.fontsize': 10, 'legend.framealpha': 0.8, 'mathtext.default': 'regular', 'axes.linewidth': 1.5, 'xtick.major.width': 1.5, 'ytick.major.width': 1.5, 'xtick.major.size': 5, 'ytick.major.size': 5, }) def identify_atom_types(struct): """原子类型识别函数""" # 初始化数据结构 atom_types = { "phosphate_oxygens": {"P-O/P=O": [], "P-OH": []}, "phosphate_hydrogens": [], "water_oxygens": [], "water_hydrogens": [], "hydronium_oxygens": [], "hydronium_hydrogens": [], "fluoride_atoms": [i for i, site in enumerate(struct) if site.species_string == "F"], "aluminum_atoms": [i for i, site in enumerate(struct) if site.species_string == "Al"] } # 构建全局KDTree all_coords = np.array([site.coords for site in struct]) kdtree = cKDTree(all_coords, boxsize=struct.lattice.abc) # 识别磷酸基团 p_atoms = [i for i, site in enumerate(struct) if site.species_string == "P"] phosphate_oxygens = [] for p_idx in p_atoms: # 查找P周围的O原子 (距离 < 1.6Å) neighbors = kdtree.query_ball_point(all_coords[p_idx], r=1.6) p_o_indices = [idx for idx in neighbors if idx != p_idx and struct[idx].species_string == "O"] phosphate_oxygens.extend(p_o_indices) # 识别所有H原子并确定归属 hydrogen_owners = {} h_atoms = [i for i, site in enumerate(struct) if site.species_string == "H"] for h_idx in h_atoms: neighbors = kdtree.query_ball_point(all_coords[h_idx], r=1.2) candidate_os = [idx for idx in neighbors if idx != h_idx and struct[idx].species_string == "O"] if not candidate_os: continue min_dist = float('inf') owner_o = None for o_idx in candidate_os: dist = struct.get_distance(h_idx, o_idx) if dist < min_dist: min_dist = dist owner_o = o_idx hydrogen_owners[h_idx] = owner_o # 分类磷酸氧:带H的为P-OH,不带H的为P-O/P=O for o_idx in phosphate_oxygens: has_hydrogen = any(owner_o == o_idx for h_idx, owner_o in hydrogen_owners.items()) if has_hydrogen: atom_types["phosphate_oxygens"]["P-OH"].append(o_idx) else: atom_types["phosphate_oxygens"]["P-O/P=O"].append(o_idx) # 识别水和水合氢离子 all_o_indices = [i for i, site in enumerate(struct) if site.species_string == "O"] non_phosphate_os = [o_idx for o_idx in all_o_indices if o_idx not in phosphate_oxygens] o_h_count = {} for h_idx, owner_o in hydrogen_owners.items(): o_h_count[owner_o] = o_h_count.get(owner_o, 0) + 1 for o_idx in non_phosphate_os: h_count = o_h_count.get(o_idx, 0) attached_hs = [h_idx for h_idx, owner_o in hydrogen_owners.items() if owner_o == o_idx] if h_count == 2: atom_types["water_oxygens"].append(o_idx) atom_types["water_hydrogens"].extend(attached_hs) elif h_count == 3: atom_types["hydronium_oxygens"].append(o_idx) atom_types["hydronium_hydrogens"].extend(attached_hs) # 识别磷酸基团的H原子 for o_idx in atom_types["phosphate_oxygens"]["P-OH"]: attached_hs = [h_idx for h_idx, owner_o in hydrogen_owners.items() if owner_o == o_idx] atom_types["phosphate_hydrogens"].extend(attached_hs) return atom_types def get_hbond_angle_config(): """氢键分析配置""" return [ { "name": "P-O/P=O···Hw", "donor_type": "water_oxygens", "acceptor_type": "P-O/P=O", "h_type": "water_hydrogens", "threshold": 2.375, "color": "blue" }, { "name": "P-O/P=O···Hh", "donor_type": "hydronium_oxygens", "acceptor_type": "P-O/P=O", "h_type": "hydronium_hydrogens", "threshold": 2.275, "color": "green" }, { "name": "P-O/P=O···Hp", "donor_type": "P-OH", "acceptor_type": "P-O/P=O", "h_type": "phosphate_hydrogens", "threshold": 2.175, "color": "red" }, { "name": "P-OH···Ow", "donor_type": "P-OH", "acceptor_type": "water_oxygens", "h_type": "phosphate_hydrogens", "threshold": 2.275, "color": "orange" }, { "name": "Hw···Ow", "donor_type": "water_oxygens", "acceptor_type": "water_oxygens", "h_type": "water_hydrogens", "threshold": 2.375, "color": "purple" }, { "name": "Hh···Ow", "donor_type": "hydronium_oxygens", "acceptor_type": "water_oxygens", "h_type": "hydronium_hydrogens", "threshold": 2.225, "color": "cyan" }, { "name": "Hw···F", "donor_type": "water_oxygens", "acceptor_type": "fluoride_atoms", "h_type": "water_hydrogens", "threshold": 2.225, "color": "magenta" }, { "name": "Hh···F", "donor_type": "hydronium_oxygens", "acceptor_type": "fluoride_atoms", "h_type": "hydronium_hydrogens", "threshold": 2.175, "color": "brown" }, { "name": "Hp···F", "donor_type": "P-OH", "acceptor_type": "fluoride_atoms", "h_type": "phosphate_hydrogens", "threshold": 2.275, "color": "pink" } ] def calculate_angle(struct, donor_idx, h_idx, acceptor_idx): """计算D-H···A键角 (角度制),使用笛卡尔向量表示并处理周期性""" # 获取分数坐标 frac_coords = struct.frac_coords lattice = struct.lattice # 获取氢原子H的分数坐标 h_frac = frac_coords[h_idx] # 计算供体D相对于H的分数坐标差 d_frac = frac_coords[donor_idx] dh_frac = d_frac - h_frac # 计算受体A相对于H的分数坐标差 a_frac = frac_coords[acceptor_idx] ah_frac = a_frac - h_frac # 应用周期性修正 (将分数坐标差限制在[-0.5, 0.5]范围内) dh_frac = np.where(dh_frac > 0.5, dh_frac - 1, dh_frac) dh_frac = np.where(dh_frac < -0.5, dh_frac + 1, dh_frac) ah_frac = np.where(ah_frac > 0.5, ah_frac - 1, ah_frac) ah_frac = np.where(ah_frac < -0.5, ah_frac + 1, ah_frac) # 转换为笛卡尔向量 (H->D 和 H->A) vec_hd = np.dot(dh_frac, lattice.matrix) # H->D向量 vec_ha = np.dot(ah_frac, lattice.matrix) # H->A向量 # 计算向量点积 dot_product = np.dot(vec_hd, vec_ha) # 计算向量模长 norm_hd = np.linalg.norm(vec_hd) norm_ha = np.linalg.norm(vec_ha) # 避免除以零 if norm_hd < 1e-6 or norm_ha < 1e-6: return None # 计算余弦值 cos_theta = dot_product / (norm_hd * norm_ha) # 处理数值误差 cos_theta = np.clip(cos_theta, -1.0, 1.0) # 计算角度 (弧度转角度) angle_rad = np.arccos(cos_theta) angle_deg = np.degrees(angle_rad) # 检查异常角度值 if angle_deg < 90 or angle_deg > 180: print(f"警告: 异常键角值 {angle_deg:.2f}° (原子 {donor_idx}-{h_idx}-{acceptor_idx})") print(f"向量HD: {vec_hd}, 长度: {norm_hd:.4f} Å") print(f"向量HA: {vec_ha}, 长度: {norm_ha:.4f} Å") print(f"点积: {dot_product:.4f}, cosθ: {cos_theta:.4f}") return angle_deg def calculate_hbond_angles_frame(struct, atom_types, hbond_config, bond_threshold=1.3): """计算单帧结构中氢键键角""" # 构建全局KDTree用于快速搜索 all_coords = np.array([site.coords for site in struct]) lattice_abc = struct.lattice.abc kdtree = cKDTree(all_coords, boxsize=lattice_abc) # 结果字典: {氢键类型: [角度列表]} angle_results = {hbond["name"]: [] for hbond in hbond_config} # 处理每一类氢键 for hbond in hbond_config: # 获取供体原子列表 if hbond["donor_type"] == "P-OH": donors = atom_types["phosphate_oxygens"]["P-OH"] else: donors = atom_types[hbond["donor_type"]] # 获取受体原子列表 if hbond["acceptor_type"] == "P-O/P=O": acceptors = atom_types["phosphate_oxygens"]["P-O/P=O"] elif hbond["acceptor_type"] == "P-OH": acceptors = atom_types["phosphate_oxygens"]["P-OH"] else: acceptors = atom_types[hbond["acceptor_type"]] # 获取氢原子列表 hydrogens = atom_types[hbond["h_type"]] # 如果没有氢原子或受体,跳过 if len(hydrogens) == 0 or len(acceptors) == 0: continue # 为受体构建KDTree(使用全局坐标) acceptor_coords = all_coords[acceptors] acceptor_kdtree = cKDTree(acceptor_coords, boxsize=lattice_abc) # 遍历所有氢原子 for h_idx in hydrogens: h_coords = all_coords[h_idx] # 查找与H成键的供体 (距离 < bond_threshold) donor_neighbors = kdtree.query_ball_point(h_coords, r=bond_threshold) donor_candidates = [idx for idx in donor_neighbors if idx in donors] # 如果没有找到供体,跳过 if not donor_candidates: continue # 选择最近的供体 min_dist = float('inf') donor_idx = None for d_idx in donor_candidates: dist = struct.get_distance(h_idx, d_idx) if dist < min_dist: min_dist = dist donor_idx = d_idx # 查找在阈值内的受体 acceptor_indices = acceptor_kdtree.query_ball_point(h_coords, r=hbond["threshold"]) for a_idx_offset in acceptor_indices: a_idx = acceptors[a_idx_offset] # 排除供体自身 if a_idx == donor_idx: continue # 计算键角 angle = calculate_angle(struct, donor_idx, h_idx, a_idx) if angle is not None and angle >= 0.0: angle_results[hbond["name"]].append(angle) return angle_results def calculate_hbond_angles_frame_wrapper(args): """包装函数用于多进程处理""" struct, hbond_config = args atom_types = identify_atom_types(struct) return calculate_hbond_angles_frame(struct, atom_types, hbond_config) def calculate_hbond_angles_parallel(structures, workers=1): """并行计算氢键键角""" hbond_config = get_hbond_angle_config() all_results = {hbond["name"]: [] for hbond in hbond_config} # 准备参数列表 args_list = [(struct, hbond_config) for struct in structures] # 如果没有可用的worker,则顺序执行 if workers == 1: results = [] for args in tqdm(args_list, desc="Calculating HBond Angles"): results.append(calculate_hbond_angles_frame_wrapper(args)) else: with multiprocessing.Pool(processes=workers) as pool: results = list(tqdm( pool.imap(calculate_hbond_angles_frame_wrapper, args_list), total=len(structures), desc="Calculating HBond Angles" )) # 合并结果 for frame_result in results: for hbond_name, angles in frame_result.items(): all_results[hbond_name].extend(angles) return all_results def plot_hbond_angle_distribution(all_angles, system_name): """绘制氢键键角分布图并保存原始数据""" # 创建输出目录 os.makedirs("HBond_Angle_Plots", exist_ok=True) os.makedirs("HBond_Angle_Data", exist_ok=True) hbond_config = get_hbond_angle_config() # 创建统计信息汇总文件 summary_path = os.path.join("HBond_Angle_Data", f"{system_name}_summary.csv") with open(summary_path, 'w', newline='') as summary_file: summary_writer = csv.writer(summary_file) summary_writer.writerow(["HBond Type", "Count", "Mean (°)", "Std Dev (°)", "Median (°)", "90% Threshold (°)", "Min (°)", "Max (°)"]) # 专业颜色方案 (JCP风格) jcp_colors = { "blue": "#1f77b4", # 深蓝色 "green": "#2ca02c", # 绿色 "red": "#d62728", # 红色 "orange": "#ff7f0e", # 橙色 "purple": "#9467bd", # 紫色 "cyan": "#17becf", # 青色 "magenta": "#e377c2", # 品红 "brown": "#8c564b", # 棕色 "pink": "#f7b6d2" # 粉红 } # 处理每种氢键类型 for hbond in hbond_config: angles = np.array(all_angles[hbond["name"]]) if len(angles) == 0: print(f"No angles found for {hbond['name']} in {system_name}") continue # 计算统计量 mean_angle = np.mean(angles) std_angle = np.std(angles) median_angle = np.median(angles) min_angle = np.min(angles) max_angle = np.max(angles) # 计算90%阈值(第10百分位数) if len(angles) >= 5: threshold_90 = np.percentile(angles, 10) else: threshold_90 = np.nan # 创建安全文件名 safe_name = hbond["name"].replace("/", "").replace("=", "").replace("···", "_").replace(" ", "_") # 保存原始角度数据 data_path = os.path.join("HBond_Angle_Data", f"{system_name}_{safe_name}.csv") with open(data_path, 'w', newline='') as data_file: data_writer = csv.writer(data_file) data_writer.writerow(["Angle (degrees)"]) for angle in angles: data_writer.writerow([f"{angle:.4f}"]) # 保存统计信息到汇总文件 with open(summary_path, 'a', newline='') as summary_file: summary_writer = csv.writer(summary_file) summary_writer.writerow([ hbond["name"], len(angles), f"{mean_angle:.4f}", f"{std_angle:.4f}", f"{median_angle:.4f}", f"{threshold_90:.4f}", f"{min_angle:.4f}", f"{max_angle:.4f}" ]) # ========== 优化绘图 ========== plt.figure(figsize=(8, 6)) # 使用JCP风格的颜色 color = jcp_colors.get(hbond["color"], hbond["color"]) # 绘制直方图 hist, bins, patches = plt.hist( angles, bins=30, density=True, alpha=0.7, color=color, edgecolor='white', linewidth=1.0, zorder=2 ) # 绘制KDE曲线 kde = gaussian_kde(angles) x = np.linspace(80, 180, 500) plt.plot( x, kde(x), color='k', linewidth=2.5, linestyle='-', zorder=3, label='KDE' ) # 添加90%阈值线 if not np.isnan(threshold_90): plt.axvline( x=threshold_90, color='#d62728', linestyle='--', linewidth=2.0, dashes=(5, 3), zorder=4, label=f'90% Threshold: {threshold_90:.1f}°' ) # ========== 优化统计信息框位置 ========== stats_text = ( f"Mean: {mean_angle:.1f}°\n" f"Median: {median_angle:.1f}°\n" f"Std Dev: {std_angle:.1f}°\n" f"Count: {len(angles)}" ) # 将统计框放在左上角图例下方 plt.text( 0.03, 0.70, # 位置:距离左边3%,距离顶部70%(图例下方) stats_text, transform=plt.gca().transAxes, fontsize=12, verticalalignment='top', horizontalalignment='left', bbox=dict( boxstyle='round', facecolor='white', alpha=0.8, edgecolor='gray', linewidth=0.8 ) ) # ========== 优化图例 ========== plt.legend( loc='upper left', # 左上角 bbox_to_anchor=(0.03, 0.98), # 距离左边3%,距离顶部2% frameon=True, framealpha=0.9, edgecolor='#333333', fancybox=False, fontsize=11 ) # 设置标题和标签 plt.title(f"{system_name}: {hbond['name']}", fontsize=16, pad=15) plt.xlabel("Bond Angle (degrees)", fontsize=14, labelpad=10) plt.ylabel("Probability Density", fontsize=14, labelpad=10) # 优化坐标轴范围 plt.xlim(80, 180) plt.ylim(0, kde(x).max() * 1.2) # 优化网格和刻度 plt.grid(True, linestyle='--', alpha=0.6, zorder=1) plt.tick_params(axis='both', which='major', labelsize=12) # 添加轻量级边框 for spine in plt.gca().spines.values(): spine.set_edgecolor('#333333') spine.set_linewidth(1.2) # 保存图像 image_path = os.path.join("HBond_Angle_Plots", f"{system_name}_{safe_name}_angle.tiff") plt.tight_layout(pad=2.0) plt.savefig(image_path, dpi=600, bbox_inches='tight') plt.close() print(f"Saved HBond angle data: {data_path}") print(f"Saved HBond angle plot: {image_path}") print(f"Saved summary statistics: {summary_path}") def main(vasprun_files, workers=1): """主处理函数""" for system_name, vasprun_file in vasprun_files.items(): print(f"\n{'='*50}") print(f"Processing {system_name}: {vasprun_file} with {workers} workers") print(f"{'='*50}") start_time = time.time() try: # 加载VASP结果 vr = Vasprun(vasprun_file, ionic_step_skip=5) structures = vr.structures print(f"Loaded {len(structures)} frames") print(f"Lattice parameters: {structures[0].lattice.abc}") print(f"Periodic boundary handling: Fractional coordinates + PBC correction") # 计算氢键键角分布 print("Calculating hydrogen bond angles...") hbond_angles = calculate_hbond_angles_parallel(structures, workers=workers) # 绘制并保存结果 plot_hbond_angle_distribution(hbond_angles, system_name) elapsed = time.time() - start_time print(f"\nCompleted processing for {system_name} in {elapsed:.2f} seconds") except Exception as e: print(f"Error processing {system_name}: {str(e)}", file=sys.stderr) import traceback traceback.print_exc() print("\nAll HBond angle analysis completed successfully!") if __name__ == "__main__": # 设置命令行参数 parser = argparse.ArgumentParser(description='Calculate hydrogen bond angles from VASP simulations') parser.add_argument('--workers', type=int, default=multiprocessing.cpu_count(), help=f'Number of parallel workers (default: {multiprocessing.cpu_count()})') args = parser.parse_args() # 自动设置vasprun文件和系统名称 vasprun_files = { "System1": "vasprun1.xml", "System2": "vasprun2.xml", "System3": "vasprun3.xml", "System4": "vasprun4.xml" } # 检查文件是否存在 missing_files = [name for name, path in vasprun_files.items() if not os.path.exists(path)] if missing_files: raise FileNotFoundError(f"Missing vasprun files: {', '.join(missing_files)}") print(f"Starting HBond angle analysis with {args.workers} workers...") main(vasprun_files, workers=args.workers)以上代码提供了氢键角度的识别以及统计,里面涉及的逻辑包括了D-H化学键的识别判定,H-A的识别判定,以及D-H-A的计算逻辑,首先必须强调的是在计算角度的时候延用笛卡尔向量坐标,角度DHA的计算由HD和HA向量计算其夹角,具体的阈值角度和距离如下def get_hbond_config(): """返回氢键配置,包含距离和角度阈值""" return [ { "name": "P-O/P=O···Hw", "donor_type": "water_oxygens", "acceptor_type": "P-O/P=O", "h_type": "water_hydrogens", "distance_threshold": 2.375, "angle_threshold": 143.99, "color": "#1f77b4" }, { "name": "P-O/P=O···Hh", "donor_type": "hydronium_oxygens", "acceptor_type": "P-O/P=O", "h_type": "hydronium_hydrogens", "distance_threshold": 2.275, "angle_threshold": 157.82, "color": "#ff7f0e" }, { "name": "P-O/P=O···Hp", "donor_type": "P-OH", "acceptor_type": "P-O/P=O", "h_type": "phosphate_hydrogens", "distance_threshold": 2.175, "angle_threshold": 155.00, "color": "#2ca02c" }, { "name": "P-OH···Ow", "donor_type": "P-OH", "acceptor_type": "water_oxygens", "h_type": "phosphate_hydrogens", "distance_threshold": 2.275, "angle_threshold": 155.13, "color": "#d62728" }, { "name": "Hw···Ow", "donor_type": "water_oxygens", "acceptor_type": "water_oxygens", "h_type": "water_hydrogens", "distance_threshold": 2.375, "angle_threshold": 138.73, "color": "#9467bd" }, { "name": "Hh···Ow", "donor_type": "hydronium_oxygens", "acceptor_type": "water_oxygens", "h_type": "hydronium_hydrogens", "distance_threshold": 2.225, "angle_threshold": 155.31, "color": "#8c564b" }, { "name": "Hw···F", "donor_type": "water_oxygens", "acceptor_type": "fluoride_atoms", "h_type": "water_hydrogens", "distance_threshold": 2.225, "angle_threshold": 137.68, "color": "#e377c2" }, { "name": "Hh···F", "donor_type": "hydronium_oxygens", "acceptor_type": "fluoride_atoms", "h_type": "hydronium_hydrogens", "distance_threshold": 2.175, "angle_threshold": 154.64, "color": "#7f7f7f" }, { "name": "Hp···F", "donor_type": "P-OH", "acceptor_type": "fluoride_atoms", "h_type": "phosphate_hydrogens", "distance_threshold": 2.275, "angle_threshold": 153.71, "color": "#bcbd22" } ] 在这里我们需要延用一样的逻辑计算不同的内容,即计算氢键寿命,采用存活相关函数SCF计算,弛豫时间采用积分法,计算对象为LAMMPS数据,并只计算一个体系,导入的文件为lammpstrij和data文件,计算的步长为0.1fs,按照合适的要求,最好符合The Journal of Chemical Physics期刊要求,轨迹输出频率,相关函数计算窗口,以及时间原点间隔,由于LAMMPS体系较大,尽可能的提高CPU和内存利用,以计算速率为主,可以增大CPU和内存的占用率。其中对原子进行映射,具体主要为1: “H”, 2: “O”, 3: “F”, 7: “P”,H3PO4-23pure.data为data文件名,nvt-P2O5-353K-23.lammpstrj为lammpstrij的文件名,然后请输出完整代码,然后提供一个测试代码,对第一帧进行氢键的识别,确保无误之后再进行氢键寿命计算的调试,在计算过程中只需要计算氢键寿命,保留键长随时间变化的识别逻辑,不需要再算键长随时间的变化。data文件为 序号 原子类型 x y z类型,例如 Atoms # charge 1 7 17.470000000000 55.085000000000 14.227000000000 2 2 16.654000000000 54.249000000000 15.498000000000 3 2 16.582000000000 55.222000000000 12.750000000000 4 2 17.569000000000 56.681000000000 14.791000000000 5 1 16.331000000000 53.417000000000 15.086000000000 优化原子ID识别问题,最好能分步进行,命令中能够输入只算第一帧氢键识别,如果无误再接着后续默认参数计算氢键寿命,然后文件的地址都设为默认参数而不用在命令栏中输入文件名,在完全实现氢键寿命的计算之前,先按照有效逻辑整理出一个能够识别氢键的代码(根据以上提供的信息,距离和角度的判别),作为测试代码,我们在这里先初步完整输出一个能够对第一帧识别所以氢键类别数量的测试代码,测试代码无误后我再要求你输出完整的关于氢键寿命的代码
最新发布
08-09
import MDAnalysis as mda import numpy as np from MDAnalysis.lib.distances import capped_distance, calc_angles, distance_array import multiprocessing as mp from tqdm import tqdm import os import argparse import matplotlib.pyplot as plt from scipy.integrate import cumtrapz import pandas as pd from functools import lru_cache import time ====================== 优化参数设置 BOX_SIZE = 80.0 CALC_REGION = [0, 60, 0, 60, 0, 60] ATOM_TYPE_MAP = {1: ‘H’, 2: ‘O’, 7: ‘P’} HB_DIST_CUTOFF = 3.5 HB_ANGLE_CUTOFF = 130 FRAME_DT = 0.1 TAU_MAX = 20.0 # 缩短时间窗口 DT_ORIGIN = 2.0 # 增加原点间隔 N_WORKERS = os.cpu_count() # 自动使用所有核心 预定义原子类型常量 O_WATER = 0 O_HYDRONIUM = 1 O_PHYDROXYL = 2 O_PDOUBLE = 3 O_PBRIDGE = 4 H_WATER = 5 H_HYDRONIUM = 6 H_PHOSPHORIC = 7 P = 8 UNKNOWN = 9 原子类型映射 ATOM_TYPE_IDS = { ‘O_water’: O_WATER, ‘O_hydronium’: O_HYDRONIUM, ‘O_phydroxyl’: O_PHYDROXYL, ‘O_pdouble’: O_PDOUBLE, ‘O_pbridge’: O_PBRIDGE, ‘H_water’: H_WATER, ‘H_hydronium’: H_HYDRONIUM, ‘H_phosphoric’: H_PHOSPHORIC, ‘P’: P, ‘UNK’: UNKNOWN } 原子类型反向映射 ATOM_TYPE_NAMES = def classify_atoms_optimized(positions, types, box): “”“向量化优化的原子分类函数”“” n_atoms = len(types) atom_labels = np.full(n_atoms, UNKNOWN, dtype=np.int8) 识别原子类型 o_indices = np.where(types == 2)[0] h_indices = np.where(types == 1)[0] p_indices = np.where(types == 7)[0] # 计算所有O-H距离 oh_distances = distance_array(positions[o_indices], positions[h_indices], box=box) oh_bonds = oh_distances < 1.2 # 计算所有P-O距离 po_distances = distance_array(positions[p_indices], positions[o_indices], box=box) po_bonds = po_distances < 1.8 # 分类氢原子 for h_idx in h_indices: # 找到与氢原子相连的氧原子 o_conn = o_indices[np.where(oh_bonds[:, np.where(h_indices == h_idx)[0][0]])[0]] if len(o_conn) == 1: o_type = atom_labels[o_conn[0]] if o_type == O_HYDRONIUM: atom_labels[h_idx] = H_HYDRONIUM elif o_type == O_WATER: atom_labels[h_idx] = H_WATER elif o_type == O_PHYDROXYL: atom_labels[h_idx] = H_PHOSPHORIC # 分类氧原子 for i, o_idx in enumerate(o_indices): # 计算连接数 h_count = np.sum(oh_bonds[i]) p_count = np.sum(po_bonds[:, i]) if h_count == 3: atom_labels[o_idx] = O_HYDRONIUM elif h_count == 2: atom_labels[o_idx] = O_WATER elif h_count == 1 and p_count == 1: atom_labels[o_idx] = O_PHYDROXYL elif p_count == 1: atom_labels[o_idx] = O_PDOUBLE elif p_count == 2: atom_labels[o_idx] = O_PBRIDGE # 分类磷原子 atom_labels[p_indices] = P return atom_labels def find_hbonds_optimized(positions, atom_labels, box): “”“向量化优化的氢键识别函数”“” 识别供体和受体 donor_mask = np.isin(atom_labels, [O_WATER, O_HYDRONIUM, O_PHYDROXYL]) acceptor_mask = np.isin(atom_labels, [O_WATER, O_HYDRONIUM, O_PDOUBLE, O_PBRIDGE, O_PHYDROXYL]) 获取供体氢原子 donor_hydrogens = [] donor_indices = [] donor_types = [] for i in np.where(donor_mask)[0]: # 找到与供体相连的氢原子 h_mask = (atom_labels == H_WATER) | (atom_labels == H_HYDRONIUM) | (atom_labels == H_PHOSPHORIC) h_positions = positions[h_mask] h_indices = np.where(h_mask)[0] dists = distance_array(positions[i].reshape(1, -1), h_positions, box=box) connected_h = np.where(dists[0] < 1.2)[0] for h_idx in connected_h: donor_hydrogens.append(h_indices[h_idx]) donor_indices.append(i) donor_types.append(atom_labels[i]) if not donor_hydrogens or not np.any(acceptor_mask): return [] # 转换为数组以提高性能 donor_hydrogens = np.array(donor_hydrogens) donor_indices = np.array(donor_indices) donor_types = np.array(donor_types) acceptor_indices = np.where(acceptor_mask)[0] # 计算供体氢到受体的距离 dists = distance_array(positions[donor_hydrogens], positions[acceptor_indices], box=box) # 找出满足距离条件的配对 valid_pairs = np.where(dists < HB_DIST_CUTOFF) hbonds = [] for i, j in zip(*valid_pairs): donor_idx = donor_indices[i] hydrogen_idx = donor_hydrogens[i] acceptor_idx = acceptor_indices[j] # 跳过自相互作用 if donor_idx == acceptor_idx: continue # 计算角度 angle = calc_angles(positions[donor_idx], positions[hydrogen_idx], positions[acceptor_idx], box=box) angle_deg = np.degrees(angle) if angle_deg > HB_ANGLE_CUTOFF: hbond_type = f"{ATOM_TYPE_NAMES[donor_types[i]].split(‘‘)[0]}-{ATOM_TYPE_NAMES[atom_labels[acceptor_idx]].split(’’)[0]}" hbonds.append({ ‘donor’: donor_idx, ‘hydrogen’: hydrogen_idx, ‘acceptor’: acceptor_idx, ‘type’: hbond_type, ‘distance’: dists[i, j], ‘angle’: angle_deg }) return hbonds def process_frame(args): “”“并行处理单个帧”“” frame_idx, top_file, traj_file, calc_region = args try: universe = mda.Universe( top_file, traj_file, topology_format=“DATA”, format=“LAMMPSDUMP”, atom_style=“id type charge x y z” ) universe.trajectory[frame_idx] 应用区域筛选 positions = universe.atoms.positions mask = ((positions[:, 0] >= calc_region[0]) & (positions[:, 0] <= calc_region[1]) & (positions[:, 1] >= calc_region[2]) & (positions[:, 1] <= calc_region[3]) & (positions[:, 2] >= calc_region[4]) & (positions[:, 2] <= calc_region[5])) if not np.any(mask): return frame_idx, [] types = universe.atoms.types.astype(int)[mask] positions = positions[mask] # 分类原子 atom_labels = classify_atoms_optimized( positions, types, universe.dimensions ) # 识别氢键 hbonds = find_hbonds_optimized( positions, atom_labels, universe.dimensions ) return frame_idx, hbonds except Exception as e: print(f"Error processing frame {frame_idx}: {str(e)}") return frame_idx, [] def calculate_hbond_lifetime(origins, universe, atom_labels_cache): “”“优化后的氢键寿命计算函数”“” n_frames = len(universe.trajectory) tau_max_frames = int(TAU_MAX / FRAME_DT) scf = np.zeros(tau_max_frames) count = np.zeros(tau_max_frames) for origin in origins: hbonds_origin = atom_labels_cache.get(origin, []) for hbond in hbonds_origin: donor = hbond[‘donor’] hydrogen = hbond[‘hydrogen’] acceptor = hbond[‘acceptor’] # 提前获取所有帧的位置 positions_cache = {} for t in range(min(tau_max_frames, n_frames - origin)): frame_idx = origin + t if frame_idx not in positions_cache: universe.trajectory[frame_idx] positions_cache[frame_idx] = universe.atoms.positions.copy() # 检查氢键存活 hbond_exists = True for t in range(1, min(tau_max_frames, n_frames - origin)): frame_idx = origin + t pos = positions_cache[frame_idx] dist = np.linalg.norm(pos[hydrogen] - pos[acceptor]) angle = calc_angles( pos[donor], pos[hydrogen], pos[acceptor], box=universe.dimensions ) angle_deg = np.degrees(angle) if dist > HB_DIST_CUTOFF or angle_deg < HB_ANGLE_CUTOFF: hbond_exists = False if hbond_exists: scf[t] += 1 count[t] += 1 else: break return scf, count def main(top_file, traj_file, output_prefix): print(f"加载拓扑文件: {top_file}“) print(f"加载轨迹文件: {traj_file}”) 创建Universe获取帧数 universe = mda.Universe( top_file, traj_file, topology_format=“DATA”, format=“LAMMPSDUMP”, atom_style=“id type charge x y z” ) n_frames = len(universe.trajectory) print(f"开始预处理 {n_frames} 帧 (并行处理)…“) start_time = time.time() # 并行处理所有帧 with mp.Pool(processes=N_WORKERS) as pool: args = [(i, top_file, traj_file, CALC_REGION) for i in range(n_frames)] results = list(tqdm(pool.imap(process_frame, args, chunksize=10), total=n_frames)) atom_labels_cache = {} all_hbonds = [] frame_hbond_counts = {} for frame_idx, frame_hbonds in results: atom_labels_cache[frame_idx] = frame_hbonds all_hbonds.extend(frame_hbonds) for hb in frame_hbonds: hb_type = hb[‘type’] if hb_type not in frame_hbond_counts: frame_hbond_counts[hb_type] = [] frame_hbond_counts[hb_type].append(1) preprocess_time = time.time() - start_time print(f"预处理完成! 耗时: {preprocess_time:.2f} 秒”) # 计算氢键密度 avg_hbonds = len(all_hbonds) / n_frames calc_volume = ((CALC_REGION[1]-CALC_REGION[0]) * (CALC_REGION[3]-CALC_REGION[2]) * (CALC_REGION[5]-CALC_REGION[4])) * 1e-27 hbond_density = avg_hbonds / (calc_volume * 6.022e23) print(f"平均氢键数: {avg_hbonds:.2f}“) print(f"氢键数密度: {hbond_density:.4f} mol/L”) # 氢键寿命计算 print(“计算氢键寿命…”) tau_max_frames = int(TAU_MAX / FRAME_DT) dt_origin_frames = int(DT_ORIGIN / FRAME_DT) origins = range(0, n_frames - tau_max_frames, dt_origin_frames) chunk_size = max(1, len(origins) // N_WORKERS) with mp.Pool(processes=N_WORKERS) as pool: tasks = [] for i in range(0, len(origins), chunk_size): chunk = origins[i:i+chunk_size] tasks.append(pool.apply_async(calculate_hbond_lifetime, (chunk, universe, atom_labels_cache))) scf_total = np.zeros(tau_max_frames) count_total = np.zeros(tau_max_frames) for task in tqdm(tasks, total=len(tasks)): scf, count = task.get() scf_total += scf count_total += count # 结果处理 scf_avg = np.where(count_total > 0, scf_total / count_total, 0) time_axis = np.arange(0, tau_max_frames) * FRAME_DT tau_relax = np.trapz(scf_avg, time_axis) print(f"弛豫时间: {tau_relax:.2f} ps") # 输出结果 scf_df = pd.DataFrame({‘Time (ps)’: time_axis, ‘SCF’: scf_avg}) scf_df.to_csv(f"{output_prefix}scf.csv", index=False) stats = { ‘Avg_HBonds_per_frame’: [avg_hbonds], ‘HBond_Density_mol/L’: [hbond_density], ‘Relaxation_Time_ps’: [tau_relax], ‘Preprocessing_Time_s’: [preprocess_time] } for hb_type, counts in frame_hbond_counts.items(): stats[f"Avg"] = [np.mean(counts)] stats_df = pd.DataFrame(stats) stats_df.to_csv(f"{output_prefix}_stats.csv", index=False) # 绘图 plt.figure(figsize=(8, 6), dpi=300) plt.plot(time_axis, scf_avg, ‘b-’, linewidth=2) plt.xlabel(‘Time (ps)’, fontsize=12, fontname=‘Times New Roman’) plt.ylabel(‘Survival Correlation Function’, fontsize=12, fontname=‘Times New Roman’) plt.title(‘Hydrogen Bond Lifetime’, fontsize=14, fontname=‘Times New Roman’) plt.grid(alpha=0.3) plt.xlim(0, TAU_MAX) plt.ylim(0, 1) plt.rcParams[‘font.family’] = ‘Times New Roman’ plt.rcParams[‘mathtext.fontset’] = ‘stix’ plt.tight_layout() plt.savefig(f"{output_prefix}_lifetime.png", format=‘png’, bbox_inches=‘tight’) print(“分析完成!”) if name == “main”: parser = argparse.ArgumentParser(description=‘优化版LAMMPS氢键寿命分析’) parser.add_argument(‘top_file’, type=str, help=‘LAMMPS拓扑文件路径 (.data)’) parser.add_argument(‘traj_file’, type=str, help=‘LAMMPS轨迹文件路径 (.lammpstrj 或 .lammpsdump)’) parser.add_argument(‘output_prefix’, type=str, help=‘输出文件前缀’) parser.add_argument(‘–workers’, type=int, default=N_WORKERS, help=‘并行工作进程数’) args = parser.parse_args() N_WORKERS = args.workers main(args.top_file, args.traj_file, args.output_prefix) 以上代码意在实现湿法磷酸中磷酸和水构成的氢键网络的氢键寿命计算,其中可能存在水合氢离子,存在动态的质子转移情况。然而在执行的时候,速度实在太慢了,目前电脑有64核 126g内存,能否进一步提高内存和CPU的利用率提高速度,总原子数在23495,体系是相当大,步长0.1fs 而总步数10000000,总时长1ns,已经能够读取data和lammpstrj文件,其中注意单位埃的输入输出,字符串类型可能会报错,为有效减小计算量,只计算前100ps内的氢键寿命.在这里你不需要为了节约时间而省略细节,尽可能的详细且设置好代码,同时命令的执行形式为python hydrogen_bond_analysis.py H3PO4-23pure.data nvt-P2O5-353K-23.lammpstrj Lifetime --workers 58,其他需要的参数放置在代码文件中谢谢
07-11
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值