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
订阅专栏 解锁全文
扫一扫
5019
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
