第十三周_Scipy题目

Exercise 10.1: Least squares 

题目描述:

Generate matrix A ∈ R m×n with m > n. Also generate some vector b ∈ Rm. Now find x = argminx||Ax−b||2. Print the norm of the residual.

代码实现:

from scipy.optimize import nnls    
import numpy as np       

m = 20
n = 10   
A = np.random.rand(m, n)    
b = np.random.rand(m)*5        
x, rnorm = nnls(A, b)

print("A = ")
print(A)
print("b = ")
print(b)   
print("rnorm = ")    
print(rnorm)

运行结果:




Exercise 10.2: Optimization

题目描述:

Find the maximum of the function:f(x) = sin2(x−2)e−x2

代码实现:

from scipy.optimize import minimize_scalar
from numpy import * 

def func(x):
	result = -(sin(x-2))**2*power(e,-x**2)
	return result

res = minimize_scalar(func)

print("The maximum = ",end="")
print(-func(res.x))

运行结果:


Exercise 10.3: Pairwise distances 

题目描述:

Let X be a matrix with n rows and m columns. How can you compute the pairwise distances between every two rows?

As an example application, consider n cities, and we are given their coordinates in two columns. Now we want a nice table that tells us for each two cities, how far they are apart.
Again, make sure you make use of Scipy’s functionality instead of writing your own routine.

代码实现:

from scipy.spatial import distance
import numpy as np    

n = 4
m = 2
X = np.random.rand(n, m)  
y = distance.pdist(X)
Z = distance.squareform(y)
print("X = ")
print(X)
print("\nMatrix Z represents the pairwise distances between two rows")
print("Z = ")
print(Z)

运行结果:


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&theta;: {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
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值