reduction_indices参数设置意思

本文详细介绍了在TensorFlow中使用tf.reduce_sum函数的方法,包括reduction_indices参数的作用及不同取值的意义,帮助读者更好地掌握该函数的使用技巧。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

在用python使用TensorFlow的时候:
tf.reduce_sum函数中reduction_indices参数表示函数的处理维度。
reduction_indices参数的值默认的时候为None,默认把所有的数据求和,即结果是一维的。
reduction_indices参数的值为0的时候,是第0维对应位置相加。
reduction_indices参数的值为1的时候,是第1维对应位置相加。
以此类推。

在这里插入图片描述

import MDAnalysis as mda import numpy as np import numba as nb from MDAnalysis.lib.distances import capped_distance, calc_angles, distance_array import multiprocessing as mp from tqdm.auto import tqdm import os import argparse import matplotlib.pyplot as plt import pandas as pd import zarr import time import gc import sys from scipy.spatial import cKDTree # ====================== 优化参数设置 ====================== DEFAULT_MAX_FRAMES = 1000 # 只处理前100ps (1000帧) DEFAULT_MEM_LIMIT = 100 # GB (126GB内存中保留6GB给系统) DEFAULT_N_WORKERS = 58 # 默认使用64个核心 BOX_SIZE = 80.0 CALC_REGION = np.array([0, 80, 0, 80, 0, 80], dtype=np.float32) ATOM_TYPE_MAP = {1: 'H', 2: 'O', 7: 'P'} HB_DIST_CUTOFF = 3.5 HB_ANGLE_CUTOFF = 130 FRAME_DT = 0.1 # ps TAU_MAX = 20.0 # ps DT_ORIGIN = 2.0 # ps MEMORY_LIMIT = 100 # GB (126GB内存中保留6GB给系统) # 原子类型常量 O_WATER, O_HYDRONIUM, O_PHYDROXYL, O_PDOUBLE, O_PBRIDGE = 0, 1, 2, 3, 4 H_WATER, H_HYDRONIUM, H_PHOSPHORIC, P, UNKNOWN = 5, 6, 7, 8, 9 # 原子类型反向映射 ATOM_TYPE_NAMES = { 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', UNKNOWN: 'UNK' } # ====================== 核心优化函数 ====================== @nb.njit(nogil=True, fastmath=True, cache=True, parallel=True) def classify_atoms_vectorized(positions, types, box): """使用Numba加速的向量化原子分类函数""" 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 = np.empty((len(o_indices), len(h_indices)), dtype=np.float32) for i in nb.prange(len(o_indices)): o_idx = o_indices[i] for j in range(len(h_indices)): h_idx = h_indices[j] d = positions[o_idx] - positions[h_idx] d -= np.round(d / box) * box # PBC校正 oh_distances[i, j] = np.sqrt(np.sum(d**2)) # 预计算P-O距离矩阵 po_distances = np.empty((len(p_indices), len(o_indices)), dtype=np.float32) for i in nb.prange(len(p_indices)): p_idx = p_indices[i] for j in range(len(o_indices)): o_idx = o_indices[j] d = positions[p_idx] - positions[o_idx] d -= np.round(d / box) * box po_distances[i, j] = np.sqrt(np.sum(d**2)) # 分类氧原子 (并行) for i in nb.prange(len(o_indices)): o_idx = o_indices[i] h_count = np.sum(oh_distances[i] < 1.2) # O-H键长阈值 p_count = 0 for j in range(len(p_indices)): if po_distances[j, i] < 1.8: # P-O键长阈值 p_count += 1 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 # 分类氢原子 (并行) for j in nb.prange(len(h_indices)): h_idx = h_indices[j] min_dist = 10.0 min_o_idx = -1 for i in range(len(o_indices)): if oh_distances[i, j] < min_dist: min_dist = oh_distances[i, j] min_o_idx = o_indices[i] if min_dist < 1.2: # O-H键长阈值 o_type = atom_labels[min_o_idx] 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 p_idx in p_indices: atom_labels[p_idx] = P return atom_labels def find_hbonds_kdtree(positions, atom_labels, box, kdtree=None): """使用KD树加速的氢键检测""" # 识别供体和受体 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_indices = np.where(donor_mask)[0] acceptor_indices = np.where(acceptor_mask)[0] if len(donor_indices) == 0 or len(acceptor_indices) == 0: return [] # 创建或复用KD树 if kdtree is None: kdtree = cKDTree(positions[acceptor_indices], boxsize=box[:3]) # 找到供体氢原子 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] hbonds = [] processed_pairs = set() # 避免重复处理 # 处理每个供体 for donor_idx in donor_indices: # 找到与供体相连的氢原子 donor_pos = positions[donor_idx] d_h_dist = distance_array(donor_pos.reshape(1, -1), h_positions, box=box)[0] connected_h = np.where(d_h_dist < 1.2)[0] for h_idx in connected_h: hydrogen_idx = h_indices[h_idx] hydrogen_pos = positions[hydrogen_idx] # 使用KD树查找附近受体 neighbors = kdtree.query_ball_point(hydrogen_pos, HB_DIST_CUTOFF) for j in neighbors: acceptor_idx = acceptor_indices[j] # 跳过自相互作用 if donor_idx == acceptor_idx: continue # 检查是否已处理 pair_id = tuple(sorted((donor_idx, acceptor_idx))) if pair_id in processed_pairs: continue processed_pairs.add(pair_id) # 修正拼写错误 # 计算角度 angle = calc_angles(donor_pos, hydrogen_pos, positions[acceptor_idx], box=box) angle_deg = np.degrees(angle) if angle_deg > HB_ANGLE_CUTOFF: dist = np.linalg.norm(hydrogen_pos - positions[acceptor_idx]) hbond_type = f"{ATOM_TYPE_NAMES[atom_labels[donor_idx]].split('_')[0]}-{ATOM_TYPE_NAMES[atom_labels[acceptor_idx]].split('_')[0]}" hbonds.append(( donor_idx, hydrogen_idx, acceptor_idx, dist, angle_deg, hbond_type )) return hbonds, kdtree # ====================== 并行处理框架 ====================== def process_frame_chunk(args): """处理一批帧的优化函数""" start_idx, end_idx, top_file, traj_file = args try: universe = mda.Universe(top_file, traj_file, topology_format="DATA", format="LAMMPSDUMP", atom_style="id type charge x y z") # 预加载区域掩码 universe.trajectory[0] 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])) filtered_indices = np.where(mask)[0] types_all = universe.atoms.types.astype(int) results = [] kdtree = None # 复用KD树 for frame_idx in range(start_idx, end_idx): if frame_idx >= len(universe.trajectory): break universe.trajectory[frame_idx] positions = universe.atoms.positions # 应用区域筛选 filtered_pos = positions[filtered_indices] filtered_types = types_all[filtered_indices] # 分类原子 atom_labels = classify_atoms_vectorized( filtered_pos, filtered_types, universe.dimensions ) # 识别氢键 (复用KD树) hbonds, kdtree = find_hbonds_kdtree( filtered_pos, atom_labels, universe.dimensions, kdtree ) # 存储压缩的氢键信息 hbond_data = [] for hb in hbonds: # 存储为元组: (donor_idx, H_idx, acceptor_idx, type_code) type_code = { 'O-O': 0, 'O-P': 1, 'P-O': 2, 'P-P': 3, 'H3O-O': 4, 'H3O-P': 5 }.get(hb[5], 6) # 默认其他类型 hbond_data.append((hb[0], hb[1], hb[2], type_code)) results.append((frame_idx, np.array(hbond_data, dtype=np.int32))) return results except Exception as e: print(f"处理帧块 [{start_idx}-{end_idx}] 时出错: {str(e)}") return [] # ====================== 氢键跟踪优化 ====================== @nb.njit(nogil=True, fastmath=True, parallel=True) def track_hbonds_batch(origin_frame, positions_cache, hbond_data, tau_max_frames, box): """批量跟踪氢键寿命的Numba优化函数""" n_hbonds = len(hbond_data) survival_matrix = np.zeros((n_hbonds, tau_max_frames), dtype=np.uint8) # 为所有氢键并行处理 for i in nb.prange(n_hbonds): donor, hydrogen, acceptor, _ = hbond_data[i] # 检查初始帧中的氢键是否存在 if donor >= positions_cache[0].shape[0] or \ hydrogen >= positions_cache[0].shape[0] or \ acceptor >= positions_cache[0].shape[0]: continue # 初始距离和角度 d0 = positions_cache[0][donor] h0 = positions_cache[0][hydrogen] a0 = positions_cache[0][acceptor] # 检查初始条件 dist0 = np.sqrt(np.sum((h0 - a0)**2)) angle_rad = calc_angles(d0, h0, a0, box=box) angle0 = np.degrees(angle_rad) if dist0 > HB_DIST_CUTOFF or angle0 < HB_ANGLE_CUTOFF: continue survival_matrix[i, 0] = 1 # 跟踪后续帧 for t in range(1, tau_max_frames): if t >= len(positions_cache): break pos_t = positions_cache[t] if donor >= pos_t.shape[0] or hydrogen >= pos_t.shape[0] or acceptor >= pos_t.shape[0]: break dt = pos_t[donor] ht = pos_t[hydrogen] at = pos_t[acceptor] # 应用PBC校正 d_vec = ht - dt d_vec -= np.round(d_vec / box) * box a_vec = at - ht a_vec -= np.round(a_vec / box) * box dist = np.sqrt(np.sum((ht - at)**2)) angle_rad = calc_angles(dt, ht, at, box=box) angle = np.degrees(angle_rad) if dist <= HB_DIST_CUTOFF and angle >= HB_ANGLE_CUTOFF: survival_matrix[i, t] = 1 else: break return survival_matrix # ====================== 氢键跟踪全局函数 ====================== def process_origin(origin, all_hbonds, traj_cache, n_frames, tau_max_frames, box): """处理单个时间原点的氢键寿命跟踪""" # 获取当前原点帧的氢键 origin_hbonds = all_hbonds[origin] if origin_hbonds is None or len(origin_hbonds) == 0: return np.zeros(tau_max_frames), np.zeros(tau_max_frames) # 修正拼写错误 # 准备帧数据缓存 frame_cache = [] for t in range(tau_max_frames): frame_idx = origin + t if frame_idx < n_frames: frame_cache.append(traj_cache.get_frame(frame_idx)) else: break # 批量跟踪氢键 survival_matrix = track_hbonds_batch( origin, frame_cache, origin_hbonds, tau_max_frames, box ) # 计算SCF scf = np.sum(survival_matrix, axis=0) count = np.array([np.sum(survival_matrix[:, :t+1].any(axis=1)) for t in range(tau_max_frames)]) return scf, count # 全局包装函数 def process_origin_wrapper(args): """包装器函数,用于传递多个参数""" return process_origin(*args) # ====================== 内存管理优化 ====================== class TrajectoryCache: """高效轨迹缓存管理器""" def __init__(self, universe, max_frames, mem_limit_gb): self.universe = universe self.n_atoms = len(universe.atoms) self.frame_size = self.n_atoms * 3 * 4 # 每帧大小 (bytes) self.max_frames = min(max_frames, len(universe.trajectory)) self.mem_limit = mem_limit_gb * 1e9 # 计算内存中可以缓存的帧数 self.cache_capacity = max(1, int(self.mem_limit * 0.7 / self.frame_size)) self.cache = {} self.access_counter = {} self.counter = 0 # 创建Zarr磁盘缓存 self.store = zarr.DirectoryStore(f'temp_positions_cache.zarr') self.root = zarr.group(store=self.store, overwrite=True) self.disk_cache = self.root.zeros( 'positions', shape=(self.max_frames, self.n_atoms, 3), chunks=(1, self.n_atoms, 3), dtype='float32' ) # 预填充磁盘缓存 for frame_idx in tqdm(range(self.max_frames), desc="预缓存轨迹"): self.universe.trajectory[frame_idx] self.disk_cache[frame_idx] = universe.atoms.positions.astype(np.float32) def get_frame(self, frame_idx): """获取帧数据,使用内存+磁盘混合缓存""" if frame_idx not in self.cache: # 如果缓存已满,移除最久未使用的 if len(self.cache) >= self.cache_capacity: least_used = min(self.access_counter, key=self.access_counter.get) del self.cache[least_used] del self.access_counter[least_used] # 从磁盘加载到内存 self.cache[frame_idx] = self.disk_cache[frame_idx][:] # 更新访问计数 self.counter += 1 self.access_counter[frame_idx] = self.counter return self.cache[frame_idx] def cleanup(self): """清理缓存""" self.cache.clear() self.access_counter.clear() gc.collect() # ====================== 主程序优化 ====================== def main(top_file, traj_file, output_prefix, max_frames=DEFAULT_MAX_FRAMES, mem_limit=DEFAULT_MEM_LIMIT, workers=DEFAULT_N_WORKERS): 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 = min(len(universe.trajectory), max_frames) n_atoms = len(universe.atoms) print(f"系统信息: {n_atoms} 个原子, {n_frames} 帧 ({n_frames*FRAME_DT}ps)") # 初始化轨迹缓存 traj_cache = TrajectoryCache(universe, n_frames, mem_limit) # ===== 第一阶段: 氢键检测 ===== print("\n===== 阶段1: 氢键检测 =====") start_time = time.time() # 计算批处理大小 frame_size = n_atoms * 3 * 4 # 位置数据大小 (bytes) batch_size = max(1, int(mem_limit * 0.3 * 1e9 / (frame_size * workers))) batches = [(i, min(i+batch_size, n_frames)) for i in range(0, n_frames, batch_size)] print(f"使用 {len(batches)} 个批次, 每批 {batch_size} 帧") print(f"启动 {workers} 个工作进程...") # 并行处理帧 all_hbonds = [None] * n_frames with mp.Pool(processes=workers) as pool: args = [(start, end, top_file, traj_file) for start, end in batches] results = list(tqdm(pool.imap_unordered(process_frame_chunk, args, chunksize=1), total=len(batches), desc="氢键检测")) # 整合结果 for batch_result in results: for frame_idx, hbonds in batch_result: all_hbonds[frame_idx] = hbonds hbond_detect_time = time.time() - start_time print(f"氢键检测完成! 耗时: {hbond_detect_time:.2f} 秒") # ===== 第二阶段: 氢键寿命跟踪 ===== print("\n===== 阶段2: 氢键寿命跟踪 =====") start_time = time.time() tau_max_frames = int(TAU_MAX / FRAME_DT) dt_origin_frames = int(DT_ORIGIN / FRAME_DT) origins = list(range(0, n_frames - tau_max_frames, dt_origin_frames)) print(f"跟踪 {len(origins)} 个时间原点的氢键寿命...") print(f"时间窗口: {TAU_MAX}ps ({tau_max_frames}帧)") # 准备结果数组 scf_total = np.zeros(tau_max_frames) count_total = np.zeros(tau_max_frames) # 准备参数列表 args_list = [(origin, all_hbonds, traj_cache, n_frames, tau_max_frames, universe.dimensions) for origin in origins] # 并行处理时间原点 - 使用全局函数 process_origin_wrapper with mp.Pool(processes=workers) as pool: results = list(tqdm( pool.imap(process_origin_wrapper, args_list, chunksize=10), total=len(origins), desc="跟踪氢键" )) for scf, count in results: valid_len = min(len(scf), len(scf_total)) scf_total[:valid_len] += scf[:valid_len] count_total[:valid_len] += count[:valid_len] # 清理缓存 traj_cache.cleanup() # ===== 结果处理 ===== print("\n===== 阶段3: 结果分析 =====") # 计算平均SCF scf_avg = np.zeros_like(scf_total) valid_mask = count_total > 0 scf_avg[valid_mask] = scf_total[valid_mask] / count_total[valid_mask] # 计算弛豫时间 time_axis = np.arange(tau_max_frames) * FRAME_DT tau_relax = np.trapz(scf_avg, time_axis) # 氢键统计 hbond_counts = {} for frame_hbonds in all_hbonds: if frame_hbonds is None: continue for hb in frame_hbonds: hb_type = hb[3] hbond_counts[hb_type] = hbond_counts.get(hb_type, 0) + 1 total_hbonds = sum(hbond_counts.values()) avg_hbonds = total_hbonds / n_frames if n_frames > 0 else 0 # 计算氢键密度 calc_volume = ((CALC_REGION[1]-CALC_REGION[0]) * (CALC_REGION[3]-CALC_REGION[2]) * (CALC_REGION[5]-CALC_REGION[4])) * 1e-27 # m³ hbond_density = avg_hbonds / (calc_volume * 6.022e23) # mol/L # 保存结果 scf_df = pd.DataFrame({ 'Time (ps)': time_axis, 'SCF': scf_avg, # 修正变量名 'Count': count_total }) scf_df.to_csv(f"{output_prefix}_scf.csv", index=False) stats_data = { 'Avg_HBonds_per_frame': avg_hbonds, 'HBond_Density_mol/L': hbond_density, 'Relaxation_Time_ps': tau_relax, 'Total_Frames': n_frames, 'HBond_Detection_Time_s': hbond_detect_time, 'Lifetime_Tracking_Time_s': time.time() - start_time } # 添加按类型的氢键统计 type_names = { 0: 'O-O', 1: 'O-P', 2: 'P-O', 3: 'P-P', 4: 'H3O-O', 5: 'H3O-P', 6: 'Other' } for type_code, count in hbond_counts.items(): stats_data[f"HBond_{type_names.get(type_code, 'Unknown')}"] = count / n_frames stats_df = pd.DataFrame([stats_data]) stats_df.to_csv(f"{output_prefix}_stats.csv", index=False) # 绘图 plt.figure(figsize=(10, 7), dpi=150) # 修正参数格式 plt.plot(time_axis, scf_avg, 'b-', linewidth=2.5, label='SCF') plt.fill_between(time_axis, scf_avg, alpha=0.2, color='blue') # 添加指数拟合 try: from scipy.optimize import curve_fit def exp_decay(t, a, tau): return a * np.exp(-t / tau) popt, _ = curve_fit(exp_decay, time_axis, scf_avg, p0=[1.0, 5.0]) plt.plot(time_axis, exp_decay(time_axis, *popt), 'r--', label=f'Exp Fit: τ={popt[1]:.2f}ps') except: pass plt.xlabel('Time (ps)', fontsize=14, fontfamily='serif') plt.ylabel('Survival Probability', fontsize=14, fontfamily='serif') plt.title(f'H-Bond Lifetime in Phosphoric Acid Solution\nτ_relax = {tau_relax:.2f} ps', fontsize=16, fontfamily='serif') plt.grid(alpha=0.3, linestyle='--') plt.xlim(0, TAU_MAX) plt.ylim(0, 1.05) plt.legend(fontsize=12) plt.tight_layout() # 添加统计信息框 stats_text = (f"Avg HBonds/frame: {avg_hbonds:.1f}\n" f"HBond density: {hbond_density:.2f} mol/L\n" f"Frames analyzed: {n_frames}") plt.figtext(0.75, 0.25, stats_text, bbox=dict(facecolor='white', alpha=0.8), fontsize=10, fontfamily='monospace') plt.savefig(f"{output_prefix}_lifetime.png", dpi=300, bbox_inches='tight') total_time = time.time() - start_time print(f"\n分析完成! 总耗时: {total_time/60:.2f} 分钟") print(f"结果保存到: {output_prefix}_*.csv/png") if __name__ == "__main__": # 简化命令行参数解析 parser = argparse.ArgumentParser( description='优化版氢键寿命分析 - 简化命令行参数', formatter_class=argparse.ArgumentDefaultsHelpFormatter ) # 必需参数 parser.add_argument('top_file', type=str, help='拓扑文件路径 (.data)') parser.add_argument('traj_file', type=str, help='轨迹文件路径 (.lammpstrj)') parser.add_argument('output_prefix', type=str, help='输出文件前缀') # 可选参数(设置默认值) parser.add_argument('--workers', type=int, default=DEFAULT_N_WORKERS, help=f'并行工作进程数 (默认: {DEFAULT_N_WORKERS})') parser.add_argument('--max_frames', type=int, default=DEFAULT_MAX_FRAMES, help=f'处理的最大帧数 (默认: {DEFAULT_MAX_FRAMES})') parser.add_argument('--mem_limit', type=int, default=DEFAULT_MEM_LIMIT, help=f'内存限制 (GB) (默认: {DEFAULT_MEM_LIMIT})') args = parser.parse_args() # 打印配置信息 print("="*60) print(f"氢键寿命分析配置:") print(f" 拓扑文件: {args.top_file}") print(f" 轨迹文件: {args.traj_file}") print(f" 输出前缀: {args.output_prefix}") print(f" 工作进程: {args.workers}") print(f" 处理帧数: {args.max_frames} (前 {args.max_frames * FRAME_DT}ps)") print(f" 内存限制: {args.mem_limit}GB") print("="*60) # 检查文件是否存在 if not os.path.exists(args.top_file): print(f"错误: 拓扑文件不存在 - {args.top_file}") sys.exit(1) if not os.path.exists(args.traj_file): print(f"错误: 轨迹文件不存在 - {args.traj_file}") sys.exit(1) # 确保输出目录存在 output_dir = os.path.dirname(args.output_prefix) if output_dir and not os.path.exists(output_dir): os.makedirs(output_dir) try: main( args.top_file, args.traj_file, args.output_prefix, max_frames=args.max_frames, mem_limit=args.mem_limit, workers=args.workers ) except Exception as e: print(f"程序执行出错: {str(e)}") import traceback traceback.print_exc() sys.exit(1) 以上代码中出现============================================================ 氢键寿命分析配置: 拓扑文件: H3PO4-23pure.data 轨迹文件: nvt-P2O5-353K-23.lammpstrj 输出前缀: Lifetime 工作进程: 58 处理帧数: 1000 (前 100.0ps) 内存限制: 100GB ============================================================ 加载拓扑文件: H3PO4-23pure.data 加载轨迹文件: nvt-P2O5-353K-23.lammpstrj 系统信息: 23659 个原子, 1000 帧 (100.0ps) 预缓存轨迹: 100%|██████████████████████████████████████████████████████████████████| 1000/1000 [01:24<00:00, 11.77it/s] ===== 阶段1: 氢键检测 ===== 使用 1 个批次, 每批 1821 帧 启动 58 个工作进程... 氢键检测: 0%| | 0/1 [00:00<?, ?it/s] 处理帧块 [0-1000] 时出错: not enough values to unpack (expected 2, got 0) 氢键检测: 100%|█████████████████████████████████████████████████████████████████████████| 1/1 [03:09<00:00, 189.14s/it] 氢键检测完成! 耗时: 189.87 秒 ===== 阶段2: 氢键寿命跟踪 ===== 跟踪 40 个时间原点的氢键寿命... 时间窗口: 20.0ps (200帧) 跟踪氢键: 0%| | 0/40 [00:00<?, ?it/s]D:\ANCONDA\envs\py13058\lib\multiprocessing\reduction.py:51: UserWarning: Reader has no dt information, set to 1.0 ps cls(buf, protocol).dump(obj) 跟踪氢键: 100%|████████████████████████████████████████████████████████████████████████| 40/40 [00:03<00:00, 10.67it/s] ===== 阶段3: 结果分析 ===== 分析完成! 总耗时: 0.13 分钟 结果保存到: Lifetime_*.csv/png
最新发布
07-12
import MDAnalysis as mda import numpy as np import numba as nb from MDAnalysis.lib.distances import capped_distance, calc_angles, distance_array import multiprocessing as mp from tqdm.auto import tqdm import os import argparse import matplotlib.pyplot as plt import pandas as pd import zarr import time import gc import sys from scipy.spatial import cKDTree # ====================== 优化参数设置 ====================== # 默认值设置为你的需求 DEFAULT_MAX_FRAMES = 1000 # 只处理前100ps (1000帧) DEFAULT_MEM_LIMIT = 100 # GB (126GB内存中保留6GB给系统) DEFAULT_N_WORKERS = 58 # 默认使用64个核心 BOX_SIZE = 80.0 CALC_REGION = np.array([0, 60, 0, 60, 0, 60], dtype=np.float32) ATOM_TYPE_MAP = {1: 'H', 2: 'O', 7: 'P'} HB_DIST_CUTOFF = 3.5 HB_ANGLE_CUTOFF = 130 FRAME_DT = 0.1 # ps TAU_MAX = 20.0 # ps DT_ORIGIN = 2.0 # ps MAX_FRAMES = 1000 # 只处理前100ps (1000帧) MEMORY_LIMIT = 100 # GB (126GB内存中保留6GB给系统) N_WORKERS = min(58, os.cpu_count()) # 使用64个核心 # 原子类型常量 O_WATER, O_HYDRONIUM, O_PHYDROXYL, O_PDOUBLE, O_PBRIDGE = 0, 1, 2, 3, 4 H_WATER, H_HYDRONIUM, H_PHOSPHORIC, P, UNKNOWN = 5, 6, 7, 8, 9 # 原子类型反向映射 ATOM_TYPE_NAMES = { 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', UNKNOWN: 'UNK' } # ====================== 核心优化函数 ====================== @nb.njit(nogil=True, fastmath=True, cache=True, parallel=True) def classify_atoms_vectorized(positions, types, box): """使用Numba加速的向量化原子分类函数""" 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 = np.empty((len(o_indices), len(h_indices)), dtype=np.float32) for i in nb.prange(len(o_indices)): o_idx = o_indices[i] for j in range(len(h_indices)): h_idx = h_indices[j] d = positions[o_idx] - positions[h_idx] d -= np.round(d / box) * box # PBC校正 oh_distances[i, j] = np.sqrt(np.sum(d**2)) # 预计算P-O距离矩阵 po_distances = np.empty((len(p_indices), len(o_indices)), dtype=np.float32) for i in nb.prange(len(p_indices)): p_idx = p_indices[i] for j in range(len(o_indices)): o_idx = o_indices[j] d = positions[p_idx] - positions[o_idx] d -= np.round(d / box) * box po_distances[i, j] = np.sqrt(np.sum(d**2)) # 分类氧原子 (并行) for i in nb.prange(len(o_indices)): o_idx = o_indices[i] h_count = np.sum(oh_distances[i] < 1.2) # O-H键长阈值 p_count = 0 for j in range(len(p_indices)): if po_distances[j, i] < 1.8: # P-O键长阈值 p_count += 1 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 # 分类氢原子 (并行) for j in nb.prange(len(h_indices)): h_idx = h_indices[j] min_dist = 10.0 min_o_idx = -1 for i in range(len(o_indices)): if oh_distances[i, j] < min_dist: min_dist = oh_distances[i, j] min_o_idx = o_indices[i] if min_dist < 1.2: # O-H键长阈值 o_type = atom_labels[min_o_idx] 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 p_idx in p_indices: atom_labels[p_idx] = P return atom_labels def find_hbonds_kdtree(positions, atom_labels, box, kdtree=None): """使用KD树加速的氢键检测""" # 识别供体和受体 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_indices = np.where(donor_mask)[0] acceptor_indices = np.where(acceptor_mask)[0] if len(donor_indices) == 0 or len(acceptor_indices) == 0: return [] # 创建或复用KD树 if kdtree is None: kdtree = cKDTree(positions[acceptor_indices], boxsize=box[:3]) # 找到供体氢原子 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] hbonds = [] processed_pairs = set() # 避免重复处理 # 处理每个供体 for donor_idx in donor_indices: # 找到与供体相连的氢原子 donor_pos = positions[donor_idx] d_h_dist = distance_array(donor_pos.reshape(1, -1), h_positions, box=box)[0] connected_h = np.where(d_h_dist < 1.2)[0] for h_idx in connected_h: hydrogen_idx = h_indices[h_idx] hydrogen_pos = positions[hydrogen_idx] # 使用KD树查找附近受体 neighbors = kdtree.query_ball_point(hydrogen_pos, HB_DIST_CUTOFF) for j in neighbors: acceptor_idx = acceptor_indices[j] # 跳过自相互作用 if donor_idx == acceptor_idx: continue # 检查是否已处理 pair_id = tuple(sorted((donor_idx, acceptor_idx))) if pair_id in processed_pairs: continue processed_pairs.add(pair_id) # 计算角度 angle = calc_angles(donor_pos, hydrogen_pos, positions[acceptor_idx], box=box) angle_deg = np.degrees(angle) if angle_deg > HB_ANGLE_CUTOFF: dist = np.linalg.norm(hydrogen_pos - positions[acceptor_idx]) hbond_type = f"{ATOM_TYPE_NAMES[atom_labels[donor_idx]].split('_')[0]}-{ATOM_TYPE_NAMES[atom_labels[acceptor_idx]].split('_')[0]}" hbonds.append(( donor_idx, hydrogen_idx, acceptor_idx, dist, angle_deg, hbond_type )) return hbonds, kdtree # ====================== 并行处理框架 ====================== def process_frame_chunk(args): """处理一批帧的优化函数""" start_idx, end_idx, top_file, traj_file = args try: universe = mda.Universe(top_file, traj_file, topology_format="DATA", format="LAMMPSDUMP", atom_style="id type charge x y z") # 预加载区域掩码 universe.trajectory[0] 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])) filtered_indices = np.where(mask)[0] types_all = universe.atoms.types.astype(int) results = [] kdtree = None # 复用KD树 for frame_idx in range(start_idx, end_idx): if frame_idx >= len(universe.trajectory): break universe.trajectory[frame_idx] positions = universe.atoms.positions # 应用区域筛选 filtered_pos = positions[filtered_indices] filtered_types = types_all[filtered_indices] # 分类原子 atom_labels = classify_atoms_vectorized( filtered_pos, filtered_types, universe.dimensions ) # 识别氢键 (复用KD树) hbonds, kdtree = find_hbonds_kdtree( filtered_pos, atom_labels, universe.dimensions, kdtree ) # 存储压缩的氢键信息 hbond_data = [] for hb in hbonds: # 存储为元组: (donor_idx, H_idx, acceptor_idx, type_code) type_code = { 'O-O': 0, 'O-P': 1, 'P-O': 2, 'P-P': 3, 'H3O-O': 4, 'H3O-P': 5 }.get(hb[5], 6) # 默认其他类型 hbond_data.append((hb[0], hb[1], hb[2], type_code)) results.append((frame_idx, np.array(hbond_data, dtype=np.int32))) return results except Exception as e: print(f"处理帧块 [{start_idx}-{end_idx}] 时出错: {str(e)}") return [] # ====================== 氢键跟踪优化 ====================== @nb.njit(nogil=True, fastmath=True, parallel=True) def track_hbonds_batch(origin_frame, positions_cache, hbond_data, tau_max_frames, box): """批量跟踪氢键寿命的Numba优化函数""" n_hbonds = len(hbond_data) survival_matrix = np.zeros((n_hbonds, tau_max_frames), dtype=np.uint8) # 为所有氢键并行处理 for i in nb.prange(n_hbonds): donor, hydrogen, acceptor, _ = hbond_data[i] # 检查初始帧中的氢键是否存在 if donor >= positions_cache[0].shape[0] or \ hydrogen >= positions_cache[0].shape[0] or \ acceptor >= positions_cache[0].shape[0]: continue # 初始距离和角度 d0 = positions_cache[0][donor] h0 = positions_cache[0][hydrogen] a0 = positions_cache[0][acceptor] # 检查初始条件 dist0 = np.sqrt(np.sum((h0 - a0)**2)) angle_rad = calc_angles(d0, h0, a0, box=box) angle0 = np.degrees(angle_rad) if dist0 > HB_DIST_CUTOFF or angle0 < HB_ANGLE_CUTOFF: continue survival_matrix[i, 0] = 1 # 跟踪后续帧 for t in range(1, tau_max_frames): if t >= len(positions_cache): break pos_t = positions_cache[t] if donor >= pos_t.shape[0] or hydrogen >= pos_t.shape[0] or acceptor >= pos_t.shape[0]: break dt = pos_t[donor] ht = pos_t[hydrogen] at = pos_t[acceptor] # 应用PBC校正 d_vec = ht - dt d_vec -= np.round(d_vec / box) * box a_vec = at - ht a_vec -= np.round(a_vec / box) * box dist = np.sqrt(np.sum((ht - at)**2)) angle_rad = calc_angles(dt, ht, at, box=box) angle = np.degrees(angle_rad) if dist <= HB_DIST_CUTOFF and angle >= HB_ANGLE_CUTOFF: survival_matrix[i, t] = 1 else: break return survival_matrix # ====================== 内存管理优化 ====================== class TrajectoryCache: """高效轨迹缓存管理器""" def __init__(self, universe, max_frames, mem_limit_gb): self.universe = universe self.n_atoms = len(universe.atoms) self.frame_size = self.n_atoms * 3 * 4 # 每帧大小 (bytes) self.max_frames = min(max_frames, len(universe.trajectory)) self.mem_limit = mem_limit_gb * 1e9 # 计算内存中可以缓存的帧数 self.cache_capacity = max(1, int(self.mem_limit * 0.7 / self.frame_size)) self.cache = {} self.access_counter = {} self.counter = 0 # 创建Zarr磁盘缓存 self.store = zarr.DirectoryStore(f'temp_positions_cache.zarr') self.root = zarr.group(store=self.store, overwrite=True) self.disk_cache = self.root.zeros( 'positions', shape=(self.max_frames, self.n_atoms, 3), chunks=(1, self.n_atoms, 3), dtype='float32' ) # 预填充磁盘缓存 for frame_idx in tqdm(range(self.max_frames), desc="预缓存轨迹"): self.universe.trajectory[frame_idx] self.disk_cache[frame_idx] = universe.atoms.positions.astype(np.float32) def get_frame(self, frame_idx): """获取帧数据,使用内存+磁盘混合缓存""" if frame_idx not in self.cache: # 如果缓存已满,移除最久未使用的 if len(self.cache) >= self.cache_capacity: least_used = min(self.access_counter, key=self.access_counter.get) del self.cache[least_used] del self.access_counter[least_used] # 从磁盘加载到内存 self.cache[frame_idx] = self.disk_cache[frame_idx][:] # 更新访问计数 self.counter += 1 self.access_counter[frame_idx] = self.counter return self.cache[frame_idx] def cleanup(self): """清理缓存""" self.cache.clear() self.access_counter.clear() gc.collect() # ====================== 主程序优化 ====================== 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 = min(len(universe.trajectory), MAX_FRAMES) n_atoms = len(universe.atoms) print(f"系统信息: {n_atoms} 个原子, {n_frames} 帧 ({n_frames*FRAME_DT}ps)") # 初始化轨迹缓存 traj_cache = TrajectoryCache(universe, n_frames, MEMORY_LIMIT) # ===== 第一阶段: 氢键检测 ===== print("\n===== 阶段1: 氢键检测 =====") start_time = time.time() # 计算批处理大小 frame_size = n_atoms * 3 * 4 # 位置数据大小 (bytes) batch_size = max(1, int(MEMORY_LIMIT * 0.3 * 1e9 / (frame_size * N_WORKERS))) batches = [(i, min(i+batch_size, n_frames)) for i in range(0, n_frames, batch_size)] print(f"使用 {len(batches)} 个批次, 每批 {batch_size} 帧") print(f"启动 {N_WORKERS} 个工作进程...") # 并行处理帧 all_hbonds = [None] * n_frames with mp.Pool(processes=N_WORKERS) as pool: args = [(start, end, top_file, traj_file) for start, end in batches] results = list(tqdm(pool.imap_unordered(process_frame_chunk, args, chunksize=1), total=len(batches), desc="氢键检测")) # 整合结果 for batch_result in results: for frame_idx, hbonds in batch_result: all_hbonds[frame_idx] = hbonds hbond_detect_time = time.time() - start_time print(f"氢键检测完成! 耗时: {hbond_detect_time:.2f} 秒") # ===== 第二阶段: 氢键寿命跟踪 ===== print("\n===== 阶段2: 氢键寿命跟踪 =====") start_time = time.time() tau_max_frames = int(TAU_MAX / FRAME_DT) dt_origin_frames = int(DT_ORIGIN / FRAME_DT) origins = list(range(0, n_frames - tau_max_frames, dt_origin_frames)) print(f"跟踪 {len(origins)} 个时间原点的氢键寿命...") print(f"时间窗口: {TAU_MAX}ps ({tau_max_frames}帧)") # 准备结果数组 scf_total = np.zeros(tau_max_frames) count_total = np.zeros(tau_max_frames) # 处理每个时间原点 def process_origin(origin): # 获取当前原点帧的氢键 origin_hbonds = all_hbonds[origin] if origin_hbonds is None or len(origin_hbonds) == 0: return np.zeros(tau_max_frames), np.zeros(tau_max_frames) # 准备帧数据缓存 frame_cache = [] for t in range(tau_max_frames): frame_idx = origin + t if frame_idx < n_frames: frame_cache.append(traj_cache.get_frame(frame_idx)) else: break # 批量跟踪氢键 survival_matrix = track_hbonds_batch( origin, frame_cache, origin_hbonds, tau_max_frames, universe.dimensions ) # 计算SCF scf = np.sum(survival_matrix, axis=0) count = np.array([np.sum(survival_matrix[:, :t+1].any(axis=1)) for t in range(tau_max_frames)]) return scf, count # 并行处理时间原点 with mp.Pool(processes=N_WORKERS) as pool: results = list(tqdm(pool.imap(process_origin, origins, chunksize=10), total=len(origins), desc="跟踪氢键")) for scf, count in results: valid_len = min(len(scf), len(scf_total)) scf_total[:valid_len] += scf[:valid_len] count_total[:valid_len] += count[:valid_len] # 清理缓存 traj_cache.cleanup() # ===== 结果处理 ===== print("\n===== 阶段3: 结果分析 =====") # 计算平均SCF scf_avg = np.zeros_like(scf_total) valid_mask = count_total > 0 scf_avg[valid_mask] = scf_total[valid_mask] / count_total[valid_mask] # 计算弛豫时间 time_axis = np.arange(tau_max_frames) * FRAME_DT tau_relax = np.trapz(scf_avg, time_axis) # 氢键统计 hbond_counts = {} for frame_hbonds in all_hbonds: if frame_hbonds is None: continue for hb in frame_hbonds: hb_type = hb[3] hbond_counts[hb_type] = hbond_counts.get(hb_type, 0) + 1 total_hbonds = sum(hbond_counts.values()) avg_hbonds = total_hbonds / n_frames if n_frames > 0 else 0 # 计算氢键密度 calc_volume = ((CALC_REGION[1]-CALC_REGION[0]) * (CALC_REGION[3]-CALC_REGION[2]) * (CALC_REGION[5]-CALC_REGION[4])) * 1e-27 # m³ hbond_density = avg_hbonds / (calc_volume * 6.022e23) # mol/L # 保存结果 scf_df = pd.DataFrame({ 'Time (ps)': time_axis, 'SCF': scf_avg, 'Count': count_total }) scf_df.to_csv(f"{output_prefix}_scf.csv", index=False) stats_data = { 'Avg_HBonds_per_frame': avg_hbonds, 'HBond_Density_mol/L': hbond_density, 'Relaxation_Time_ps': tau_relax, 'Total_Frames': n_frames, 'HBond_Detection_Time_s': hbond_detect_time, 'Lifetime_Tracking_Time_s': time.time() - start_time } # 添加按类型的氢键统计 type_names = { 0: 'O-O', 1: 'O-P', 2: 'P-O', 3: 'P-P', 4: 'H3O-O', 5: 'H3O-P', 6: 'Other' } for type_code, count in hbond_counts.items(): stats_data[f"HBond_{type_names.get(type_code, 'Unknown')}"] = count / n_frames stats_df = pd.DataFrame([stats_data]) stats_df.to_csv(f"{output_prefix}_stats.csv", index=False) # 绘图 plt.figure(figsize=(10, 7), dpi=150) plt.plot(time_axis, scf_avg, 'b-', linewidth=2.5, label='SCF') plt.fill_between(time_axis, scf_avg, alpha=0.2, color='blue') # 添加指数拟合 try: from scipy.optimize import curve_fit def exp_decay(t, a, tau): return a * np.exp(-t / tau) popt, _ = curve_fit(exp_decay, time_axis, scf_avg, p0=[1.0, 5.0]) plt.plot(time_axis, exp_decay(time_axis, *popt), 'r--', label=f'Exp Fit: τ={popt[1]:.2f}ps') except: pass plt.xlabel('Time (ps)', fontsize=14, fontfamily='serif') plt.ylabel('Survival Probability', fontsize=14, fontfamily='serif') plt.title(f'H-Bond Lifetime in Phosphoric Acid Solution\nτ_relax = {tau_relax:.2f} ps', fontsize=16, fontfamily='serif') plt.grid(alpha=0.3, linestyle='--') plt.xlim(0, TAU_MAX) plt.ylim(0, 1.05) plt.legend(fontsize=12) plt.tight_layout() # 添加统计信息框 stats_text = (f"Avg HBonds/frame: {avg_hbonds:.1f}\n" f"HBond density: {hbond_density:.2f} mol/L\n" f"Frames analyzed: {n_frames}") plt.figtext(0.75, 0.25, stats_text, bbox=dict(facecolor='white', alpha=0.8), fontsize=10, fontfamily='monospace') plt.savefig(f"{output_prefix}_lifetime.png", dpi=300, bbox_inches='tight') total_time = time.time() - start_time print(f"\n分析完成! 总耗时: {total_time/60:.2f} 分钟") print(f"结果保存到: {output_prefix}_*.csv/png") if __name__ == "__main__": # 简化命令行参数解析 parser = argparse.ArgumentParser( description='优化版氢键寿命分析 - 简化命令行参数', formatter_class=argparse.ArgumentDefaultsHelpFormatter ) # 必需参数 parser.add_argument('top_file', type=str, help='拓扑文件路径 (.data)') parser.add_argument('traj_file', type=str, help='轨迹文件路径 (.lammpstrj)') parser.add_argument('output_prefix', type=str, help='输出文件前缀') # 可选参数(设置默认值) parser.add_argument('--workers', type=int, default=DEFAULT_N_WORKERS, help=f'并行工作进程数 (默认: {DEFAULT_N_WORKERS})') # 隐藏其他参数(使用默认值) parser.add_argument('--max_frames', type=int, default=DEFAULT_MAX_FRAMES, help=argparse.SUPPRESS) # 隐藏参数,但仍可使用 parser.add_argument('--mem_limit', type=int, default=DEFAULT_MEM_LIMIT, help=argparse.SUPPRESS) # 隐藏参数,但仍可使用 args = parser.parse_args() # 打印配置信息 print("="*60) print(f"氢键寿命分析配置:") print(f" 拓扑文件: {args.top_file}") print(f" 轨迹文件: {args.traj_file}") print(f" 输出前缀: {args.output_prefix}") print(f" 工作进程: {args.workers}") print(f" 处理帧数: {args.max_frames} (前 {args.max_frames * FRAME_DT}ps)") print(f" 内存限制: {args.mem_limit}GB") print("="*60) # 检查文件是否存在 if not os.path.exists(args.top_file): print(f"错误: 拓扑文件不存在 - {args.top_file}") sys.exit(1) if not os.path.exists(args.traj_file): print(f"错误: 轨迹文件不存在 - {args.traj_file}") sys.exit(1) # 确保输出目录存在 output_dir = os.path.dirname(args.output_prefix) if output_dir and not os.path.exists(output_dir): os.makedirs(output_dir) try: main(args.top_file, args.traj_file, args.output_prefix) except Exception as e: print(f"程序执行出错: {str(e)}") import traceback traceback.print_exc() sys.exit(1)这个代码目前可运行,但是报错Traceback (most recent call last): File "hydrogen_bond_analysis6.py", line 634, in <module> main(args.top_file, args.traj_file, args.output_prefix) File "hydrogen_bond_analysis6.py", line 475, in main results = list(tqdm(pool.imap(process_origin, origins, chunksize=10), File "D:\ANCONDA\envs\py13058\lib\site-packages\tqdm\std.py", line 1181, in __iter__ for obj in iterable: File "D:\ANCONDA\envs\py13058\lib\multiprocessing\pool.py", line 420, in <genexpr> return (item for chunk in result for item in chunk) File "D:\ANCONDA\envs\py13058\lib\multiprocessing\pool.py", line 868, in next raise value File "D:\ANCONDA\envs\py13058\lib\multiprocessing\pool.py", line 537, in _handle_tasks put(task) File "D:\ANCONDA\envs\py13058\lib\multiprocessing\connection.py", line 206, in send self._send_bytes(_ForkingPickler.dumps(obj)) File "D:\ANCONDA\envs\py13058\lib\multiprocessing\reduction.py", line 51, in dumps cls(buf, protocol).dump(obj) AttributeError: Can't pickle local object 'main.<locals>.process_origin'
07-11
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值