Ionic 2.0 Alpha版本发布了

Ionic团队宣布发布Ionic 2.0 Alpha版本,在Angular Connect会议上首次公开。该版本包含许多重大改进,如增强导航体验、增加原生功能支持、强大的主题定制能力等。

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

Today, at the Angular Connect conference in London, we announced the first public alpha version of Ionic 2. Many months of hard work, false starts, lessons learned, and “aha!” moments have led us to this day, and we are so excited to finally take the wraps off the next generation of Ionic.

We’ve talked a lot about why we built Ionic 2, but it’s worth revisiting the story now that Ionic 2 actually exists in the wild. We started working on Ionic in the Summer of 2013, back when iOS 6 was king and Android 2.3 was still powering a significant number of devices. We had different browser APIs (lack thereof, really), and browser engines that seriously underperform what we have today. We decided to base Ionic on Angular, which ended up being a really lucky decision, but Angular 1.1 which we used back then is a far cry from new Angular 1 versions, and a generation in capability from Angular 2.

Fast forward two years, and Ionic has found widespread adoption in a diverse set of industries from developers around the world. Collectively, over 1.2M apps have been built on Ionic, a number that is increasing exponentially. With all these developers creating all these apps on Ionic, we’ve been able to learn a thing or two about how to build the best possible mobile development toolkit to help people build great apps quickly, without breaking the bank.

With all that data, we realized there were a number of crucial things we could improve in Ionic, but they would require more significant architectural changes. Angular 2 seemed like the perfect opportunity to do just that.

With Ionic 2, we’ve overhauled and added a number of important features, including:

  • Overhauled Navigation: Completely control the navigation experience of your app without being tied to the URL bar. Navigate to any page inside of any view, including modals, side menus, and other view containers, while maintaining full deeplinking capability. This will completely change how you think about your app experience.

  • Native Support: We’ve added more native functionality directly into Ionic, making it easy to take advantage of the full power of the device without hunting down external plugins and code.

  • Powerful Theming: Don’t build apps that look like stock iOS/Android/Ionic. With the new theming system, it’s easy to instantly match your brand colors and design.

  • Material Design: Full material design support for Android apps.

Not to mention new and improved components, new docs, a new animation system, dramatically improved performance, and a whole lot more.

Let’s address the elephant in the room: “What about Ionic 1? What’s going to happen to it?” First of all, we will absolutely continue supporting Ionic 1 for a long time. In fact, there are a number of improvements we’re working on, including upgrading to the next version of Angular. At the same time, we are going to work to make it easy to upgrade to Ionic 2 slowly over time, as the community also starts to adopt Angular 2. While we aren’t going to apologize for wanting to push Ionic hard into the future, we also understand that software can live for a long time, and should be supported for as long as possible.

One of the things that gets us so excited about Ionic 2 is just how much room there was to improve on Ionic 1, especially when it came to performance and native functionality. We’ve been able to really optimize how Ionic apps feel with v2, and we think it’s going to completely change what you thought the web was capable of on mobile.

To get started, visit ionic.io/2, or visit the Getting Started page. We’ve also opened up a temporary GitHub Repository for Ionic 2 (see the README for a few demos!), and we will be moving the code over to the official Ionic repo over the coming weeks. If you run into any issues or have any feedback, please file issues on that repo! Also check out Adam’s great slides from his Ionic 2 talk today, complete with demos!

We want to thank everyone from the Ionic community for your support over the last two years. None of us expected that Ionic would become the most popular cross-platform mobile development toolkit in such a short amount of time, and it’s really your passion for building apps and sharing that passion with the world that has made any of this possible.

From all of us at Ionic, thank you, and we look forward to your feedback on Ionic 2 as we all work hard to push open web technologies on mobile into the future!

你先读取该代码import numpy as np import matplotlib.pyplot as plt from pymatgen.io.vasp import Vasprun from pymatgen.core.structure import Structure from scipy.signal import savgol_filter from scipy.spatial import cKDTree from tqdm import tqdm import matplotlib as mpl import warnings from collections import defaultdict import os import csv import argparse import multiprocessing from functools import partial import time import dill # 忽略可能的警告 warnings.filterwarnings("ignore", category=UserWarning) # 专业绘图设置 - 符合Journal of Chemical Physics要求 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, }) # 1. 增强的原子类型识别函数 - 逐帧识别 def identify_atom_types(struct): """识别所有关键原子类型并排除自身化学键""" # 磷酸氧分类 p_oxygens = {"P=O": [], "P-O": [], "P-OH": []} phosphate_hydrogens = [] # 仅P-OH基团中的H原子 # 水合氢离子识别 hydronium_oxygens = [] hydronium_hydrogens = [] # H₃O⁺中的H原子 # 普通水分子 water_oxygens = [] water_hydrogens = [] # 普通水中的H原子 # 氟离子 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"] # 创建快速邻居查找表 neighbor_cache = defaultdict(list) for i, site in enumerate(struct): if site.species_string == "O": neighbors = struct.get_neighbors(site, r=1.3) h_neighbors = [n[0] for n in neighbors if n[0].species_string == "H"] neighbor_cache[i] = h_neighbors # 识别水合氢离子 (H₃O⁺) if len(h_neighbors) == 3: hydronium_oxygens.append(i) for h_site in h_neighbors: hydronium_hydrogens.append(h_site.index) # 识别磷酸基团 for site in struct: if site.species_string == "P": neighbors = struct.get_neighbors(site, r=2.0) # 扩大搜索半径 # 筛选氧原子邻居 o_neighbors = [(n[0], n[1]) for n in neighbors if n[0].species_string == "O"] if len(o_neighbors) < 4: # 如果找不到4个氧原子,使用旧方法 for neighbor in o_neighbors: nn_site = neighbor[0] if neighbor[1] < 1.55: p_oxygens["P=O"].append(nn_site.index) else: if any(n[0].species_string == "H" for n in struct.get_neighbors(nn_site, r=1.3)): p_oxygens["P-OH"].append(nn_site.index) else: p_oxygens["P-O"].append(nn_site.index) continue # 按距离排序 o_neighbors.sort(key=lambda x: x[1]) # 最近的氧原子为P=O p_double_o = o_neighbors[0][0] p_oxygens["P=O"].append(p_double_o.index) # 其他三个氧原子 for i in range(1, 4): o_site = o_neighbors[i][0] # 检查氧原子上是否有氢 if neighbor_cache.get(o_site.index, []): p_oxygens["P-OH"].append(o_site.index) else: p_oxygens["P-O"].append(o_site.index) # 识别P-OH基团中的H原子 (磷酸中的H) for o_idx in p_oxygens["P-OH"]: # 获取与P-OH氧相连的H原子 h_neighbors = neighbor_cache.get(o_idx, []) for h_site in h_neighbors: if h_site.species_string == "H": phosphate_hydrogens.append(h_site.index) # 识别普通水分子 (排除磷酸氧和水合氢离子) for i, site in enumerate(struct): if site.species_string == "O" and i not in hydronium_oxygens: is_phosphate_oxygen = False for cat in p_oxygens.values(): if i in cat: is_phosphate_oxygen = True break if not is_phosphate_oxygen: water_oxygens.append(i) # 识别普通水分子中的H原子 (水中的H) for o_idx in water_oxygens: h_neighbors = neighbor_cache.get(o_idx, []) for h_site in h_neighbors: if h_site.species_string == "H": water_hydrogens.append(h_site.index) return { "phosphate_oxygens": p_oxygens, "phosphate_hydrogens": phosphate_hydrogens, "water_oxygens": water_oxygens, "water_hydrogens": water_hydrogens, "hydronium_oxygens": hydronium_oxygens, "hydronium_hydrogens": hydronium_hydrogens, "fluoride_atoms": fluoride_atoms, "aluminum_atoms": aluminum_atoms } # 2. RDF计算函数 - 修复负值问题和序列化问题 def process_frame(struct, center_sel, target_sel, r_max, exclude_bonds, bond_threshold): """处理单帧结构计算,完全处理空原子类型情况""" # 每帧重新识别原子类型(关键!) atom_types = identify_atom_types(struct) # 获取中心原子和目标原子 centers = center_sel(atom_types) targets = target_sel(atom_types) # 处理空原子类型情况 - 第一重保护 if len(centers) == 0 or len(targets) == 0: return { "distances": np.array([], dtype=np.float64), "n_centers": 0, "n_targets": 0, "volume": struct.volume } center_coords = np.array([struct[i].coords for i in centers]) target_coords = np.array([struct[i].coords for i in targets]) lattice = struct.lattice kdtree = cKDTree(target_coords, boxsize=lattice.abc) # 动态确定邻居数量 - 不超过目标原子数 k_val = min(50, len(targets)) # 处理目标原子数量为0的情况 - 第二重保护 if k_val == 0: return { "distances": np.array([], dtype=np.float64), "n_centers": len(centers), "n_targets": len(targets), "volume": struct.volume } # 执行查询并确保结果统一格式 try: query_result = kdtree.query(center_coords, k=k_val, distance_upper_bound=r_max) except Exception as e: # 异常处理 - 返回空结果 print(f"KDTree query error: {str(e)}") return { "distances": np.array([], dtype=np.float64), "n_centers": len(centers), "n_targets": len(targets), "volume": struct.volume } # 统一处理不同维度的返回结果 if k_val == 1: # 处理单邻居情况 if isinstance(query_result, tuple): distances, indices = query_result else: distances = query_result indices = np.zeros_like(distances, dtype=int) # 确保数组格式 distances = np.atleast_1d(distances) indices = np.atleast_1d(indices) else: # 多邻居情况 distances, indices = query_result # 确保二维数组格式 if distances.ndim == 1: distances = distances.reshape(-1, 1) indices = indices.reshape(-1, 1) valid_distances = [] for i in range(distances.shape[0]): center_idx = centers[i] for j in range(distances.shape[1]): dist = distances[i, j] # 跳过超出范围的距离 if dist > r_max or np.isinf(dist): continue target_idx = targets[indices[i, j]] # 排除化学键 if exclude_bonds: actual_dist = struct.get_distance(center_idx, target_idx) if actual_dist < bond_threshold: continue valid_distances.append(dist) return { "distances": np.array(valid_distances, dtype=np.float64), "n_centers": len(centers), "n_targets": len(targets), "volume": struct.volume } def calculate_rdf_parallel(structures, center_sel, target_sel, r_max=8.0, bin_width=0.05, exclude_bonds=True, bond_threshold=1.3, workers=1): """ 并行计算径向分布函数 :param workers: 并行工作进程数 """ bins = np.arange(0, r_max, bin_width) hist = np.zeros(len(bins) - 1) total_centers = 0 total_targets = 0 total_volume = 0 # 准备参数 - 使用dill解决序列化问题 dill.settings['recurse'] = True func = partial(process_frame, center_sel=center_sel, target_sel=target_sel, r_max=r_max, exclude_bonds=exclude_bonds, bond_threshold=bond_threshold) # 使用多进程池 with multiprocessing.Pool(processes=workers) as pool: results = [] # 使用imap_unordered提高效率 for res in tqdm(pool.imap_unordered(func, structures), total=len(structures), desc="Calculating RDF"): results.append(res) # 处理结果 - 特别注意空结果处理 n_frames = 0 for res in results: if res is None: continue n_frames += 1 valid_distances = res["distances"] n_centers = res["n_centers"] n_targets = res["n_targets"] volume = res["volume"] # 累加计数 if len(valid_distances) > 0: hist += np.histogram(valid_distances, bins=bins)[0] total_centers += n_centers total_targets += n_targets total_volume += volume # 修正归一化 - 解决负值问题 if n_frames == 0: # 没有有效帧时返回空结果 r = bins[:-1] + bin_width/2 return r, np.zeros_like(r), {"position": None, "value": None} avg_density = total_targets / total_volume if total_volume > 0 else 0 r = bins[:-1] + bin_width/2 rdf = np.zeros_like(r) for i in range(len(hist)): r_lower = bins[i] r_upper = bins[i+1] shell_vol = 4/3 * np.pi * (r_upper**3 - r_lower**3) expected_count = shell_vol * avg_density * total_centers # 避免除以零 if expected_count > 1e-10: rdf[i] = hist[i] / expected_count else: rdf[i] = 0 # 更稳健的平滑处理 - 避免边界效应 if len(rdf) > 10: window_length = min(15, len(rdf)//2*2+1) polyorder = min(5, window_length-1) rdf_smoothed = savgol_filter(rdf, window_length=window_length, polyorder=polyorder, mode='mirror') else: rdf_smoothed = rdf # 计算主要峰值 peak_info = {} mask = (r >= 1.5) & (r <= 3.0) if np.any(mask) and np.any(rdf_smoothed[mask] > 0): peak_idx = np.argmax(rdf_smoothed[mask]) peak_pos = r[mask][peak_idx] peak_val = rdf_smoothed[mask][peak_idx] peak_info = {"position": peak_pos, "value": peak_val} else: peak_info = {"position": None, "value": None} return r, rdf_smoothed, peak_info # 3. 定义精细化的选择器函数(避免lambda序列化问题) def selector_phosphate_P_double_O(atom_types): return atom_types["phosphate_oxygens"]["P=O"] def selector_phosphate_P_OH(atom_types): return atom_types["phosphate_oxygens"]["P-OH"] def selector_phosphate_P_O(atom_types): return atom_types["phosphate_oxygens"]["P-O"] def selector_phosphate_hydrogens(atom_types): return atom_types["phosphate_hydrogens"] def selector_water_only_hydrogens(atom_types): """仅选择水分子中的氢原子""" return atom_types["water_hydrogens"] def selector_hydronium_only_hydrogens(atom_types): """仅选择水合氢离子中的氢原子""" return atom_types["hydronium_hydrogens"] def selector_water_only_oxygens(atom_types): """仅选择水分子中的氧原子""" return atom_types["water_oxygens"] def selector_hydronium_only_oxygens(atom_types): """仅选择水合氢离子中的氧原子""" return atom_types["hydronium_oxygens"] def selector_fluoride_atoms(atom_types): return atom_types["fluoride_atoms"] def selector_aluminum_atoms(atom_types): return atom_types["aluminum_atoms"] def selector_all_phosphate_oxygens(atom_types): return (atom_types["phosphate_oxygens"]["P=O"] + atom_types["phosphate_oxygens"]["P-O"] + atom_types["phosphate_oxygens"]["P-OH"]) # 4. 根据您的要求定义六张图的RDF分组配置 def get_rdf_groups(): """返回六张图的RDF分组配置(完全符合您的需求)""" return { # 图1: Al的配位情况 "Al_Coordination": [ (selector_aluminum_atoms, selector_fluoride_atoms, "Al-F", "blue"), (selector_aluminum_atoms, selector_water_only_oxygens, "Al-Ow", "green"), (selector_aluminum_atoms, selector_all_phosphate_oxygens, "Al-Op", "red") ], # 图2: F与H形成的氢键 "F_Hydrogen_Bonding": [ (selector_fluoride_atoms, selector_water_only_hydrogens, "F-Hw", "lightblue"), (selector_fluoride_atoms, selector_hydronium_only_hydrogens, "F-Hh", "blue"), (selector_fluoride_atoms, selector_phosphate_hydrogens, "F-Hp", "darkblue") ], # 图3: 磷酸作为受体与周围环境的氢键(区分氧类型) "Phosphate_Acceptor": [ (selector_phosphate_P_double_O, selector_water_only_hydrogens, "P=O···Hw", "orange"), (selector_phosphate_P_double_O, selector_hydronium_only_hydrogens, "P=O···Hh", "red"), (selector_phosphate_P_O, selector_water_only_hydrogens, "P-O···Hw", "lightgreen"), (selector_phosphate_P_O, selector_hydronium_only_hydrogens, "P-O···Hh", "green"), (selector_phosphate_P_OH, selector_water_only_hydrogens, "P-OH···Hw", "lightblue"), (selector_phosphate_P_OH, selector_hydronium_only_hydrogens, "P-OH···Hh", "blue") ], # 图4: 磷酸-水-水合氢离子交叉氢键(排除同种类型) "Cross_Species_HBonding": [ (selector_phosphate_hydrogens, selector_water_only_oxygens, "Hp···Ow", "pink"), (selector_phosphate_hydrogens, selector_hydronium_only_oxygens, "Hp···Oh", "purple"), (selector_water_only_hydrogens, selector_all_phosphate_oxygens, "Hw···Op", "lightgreen"), (selector_water_only_hydrogens, selector_hydronium_only_oxygens, "Hw···Oh", "green"), (selector_hydronium_only_hydrogens, selector_water_only_oxygens, "Hh···Ow", "lightblue"), (selector_hydronium_only_hydrogens, selector_all_phosphate_oxygens, "Hh···Op", "blue") ], # 图5: 同类型分子内/间氢键(区分磷酸氧类型) "Same_Species_HBonding": [ (selector_phosphate_hydrogens, selector_phosphate_P_double_O, "Hp···P=O", "red"), (selector_phosphate_hydrogens, selector_phosphate_P_O, "Hp···P-O", "orange"), (selector_phosphate_hydrogens, selector_phosphate_P_OH, "Hp···P-OH", "yellow"), (selector_water_only_hydrogens, selector_water_only_oxygens, "Hw···Ow", "lightblue"), (selector_hydronium_only_hydrogens, selector_hydronium_only_oxygens, "Hh···Oh", "blue") ], # 图6: O-O聚集分析(Op不区分类型) "O_O_Aggregation": [ (selector_all_phosphate_oxygens, selector_water_only_oxygens, "Op-Ow", "blue"), (selector_all_phosphate_oxygens, selector_hydronium_only_oxygens, "Op-Oh", "green"), (selector_all_phosphate_oxygens, selector_all_phosphate_oxygens, "Op-Op", "red"), (selector_water_only_oxygens, selector_hydronium_only_oxygens, "Ow-Oh", "purple"), (selector_water_only_oxygens, selector_water_only_oxygens, "Ow-Ow", "cyan"), (selector_hydronium_only_oxygens, selector_hydronium_only_oxygens, "Oh-Oh", "magenta") ] } # 5. 主程序 - 优化并行处理 def main(workers=1): # 定义要处理的体系 vasprun_files = { "System1": "vasprun1.xml", "System2": "vasprun2.xml", "System3": "vasprun3.xml", "System4": "vasprun4.xml" } # 获取RDF分组配置 rdf_groups = get_rdf_groups() # 标题映射(根据您的要求) title_map = { "Al_Coordination": "Al Coordination Environment", "F_Hydrogen_Bonding": "F-H Hydrogen Bonding", "Phosphate_Acceptor": "Phosphate as H-bond Acceptor", "Cross_Species_HBonding": "Cross H-bonding between Different Species", "Same_Species_HBonding": "Intra- and Inter-molecular H-bonding", "O_O_Aggregation": "O-O Aggregation Analysis" } # 存储所有数据 all_system_data = {} group_y_max = {group_name: 0 for group_name in list(rdf_groups.keys())} group_x_max = { "Al_Coordination": (1.5, 3.5), "F_Hydrogen_Bonding": (1.0, 3.0), "Phosphate_Acceptor": (1.0, 3.0), "Cross_Species_HBonding": (1.0, 3.0), "Same_Species_HBonding": (1.0, 3.0), "O_O_Aggregation": (2.0, 6.0) } # 创建输出目录 os.makedirs("RDF_Plots", exist_ok=True) # 计算所有体系的所有RDF数据 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") # 存储体系数据 system_data = { "rdf_results": {}, "peak_infos": {} } # 计算所有RDF分组 for group_name, pairs in rdf_groups.items(): system_data["rdf_results"][group_name] = {} system_data["peak_infos"][group_name] = {} group_y_max_current = 0 for center_sel, target_sel, label, color in pairs: print(f"\nCalculating RDF for: {label}") try: r, rdf, peak_info = calculate_rdf_parallel( structures, center_sel, target_sel, r_max=10.0, exclude_bonds=True, bond_threshold=1.3, workers=workers ) system_data["rdf_results"][group_name][label] = (r, rdf, color) system_data["peak_infos"][group_name][label] = peak_info if len(rdf) > 0: current_max = np.max(rdf) if current_max > group_y_max_current: group_y_max_current = current_max if peak_info["position"] is not None: print(f" Peak for {label}: {peak_info['position']:.3f} Å (g(r) = {peak_info['value']:.2f})") else: print(f" No significant peak found for {label} in 1.5-3.0 Å range") except Exception as e: print(f"Error calculating RDF for {label}: {str(e)}") system_data["rdf_results"][group_name][label] = (np.array([]), np.array([]), color) system_data["peak_infos"][group_name][label] = {"position": None, "value": None} if group_y_max_current > group_y_max[group_name]: group_y_max[group_name] = group_y_max_current all_system_data[system_name] = system_data 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)}") # 为每个分组添加余量 for group_name in group_y_max: group_y_max[group_name] = max(group_y_max[group_name] * 1.15, 3.0) # 确保最小值 # 第二步:生成符合期刊要求的图表 for system_name, system_data in all_system_data.items(): print(f"\nGenerating publication-quality plots for {system_name}") for group_name, group_data in system_data["rdf_results"].items(): fig, ax = plt.subplots(figsize=(8, 6)) # 设置坐标轴范围 xlim = group_x_max.get(group_name, (0, 6.0)) ylim = (0, group_y_max[group_name]) for label, (r, rdf, color) in group_data.items(): if len(r) > 0 and len(rdf) > 0: ax.plot(r, rdf, color=color, label=label, linewidth=2.0) ax.set_xlim(xlim) ax.set_ylim(ylim) # 期刊格式标签 ax.set_xlabel('Radial Distance (Å)', fontweight='bold') ax.set_ylabel('g(r)', fontweight='bold') # 添加体系名称到标题 ax.set_title(f"{system_name}: {title_map[group_name]}", fontsize=16, pad=15) # 精简图例 ncol = 3 if group_name == "Same_Species_HBonding" else 1 # 图5使用三列图例 ax.legend(ncol=ncol, loc='best', framealpha=0.8, fontsize=10) # 添加氢键区域标记(除O-O聚集图外) if group_name != "O_O_Aggregation": ax.axvspan(1.5, 2.5, alpha=0.1, color='green', zorder=0) ax.text(1.7, ylim[1]*0.85, 'H-bond Region', fontsize=10) # 添加网格 ax.grid(True, linestyle='--', alpha=0.5) # 保存高分辨率图片 plt.tight_layout() filename = os.path.join("RDF_Plots", f"RDF_{system_name}_{group_name}.tiff") plt.savefig(filename, bbox_inches='tight', dpi=600, format='tiff') print(f"Saved publication plot: {filename}") plt.close() # 保存Origin兼容数据 save_origin_data(system_name, system_data) print("\nAll RDF analysis completed successfully!") def save_origin_data(system_name, system_data): """保存Origin兼容格式数据""" os.makedirs("Origin_Data", exist_ok=True) system_dir = os.path.join("Origin_Data", system_name) os.makedirs(system_dir, exist_ok=True) # 保存峰值信息 peak_info_path = os.path.join(system_dir, f"Peak_Positions_{system_name}.csv") with open(peak_info_path, 'w', newline='') as csvfile: writer = csv.writer(csvfile) writer.writerow(["Group", "Interaction", "Peak Position (A)", "g(r) Value"]) for group_name, peaks in system_data["peak_infos"].items(): for label, info in peaks.items(): if info["position"] is not None: writer.writerow([group_name, label, f"{info['position']:.3f}", f"{info['value']:.3f}"]) else: writer.writerow([group_name, label, "N/A", "N/A"]) print(f"Saved peak positions: {peak_info_path}") # 保存RDF数据 for group_name, group_results in system_data["rdf_results"].items(): group_dir = os.path.join(system_dir, group_name) os.makedirs(group_dir, exist_ok=True) for label, (r, rdf, color) in group_results.items(): if len(r) > 0 and len(rdf) > 0: safe_label = label.replace(" ", "_").replace("/", "_").replace("=", "_") safe_label = safe_label.replace("(", "").replace(")", "").replace("$", "") filename = f"RDF_{system_name}_{group_name}_{safe_label}.csv" filepath = os.path.join(group_dir, filename) with open(filepath, 'w', newline='') as csvfile: writer = csv.writer(csvfile) writer.writerow(["Distance (A)", "g(r)"]) for i in range(len(r)): writer.writerow([f"{r[i]:.6f}", f"{rdf[i]:.6f}"]) print(f"Saved Origin data: {filename}") if __name__ == "__main__": # 设置命令行参数 parser = argparse.ArgumentParser(description='Calculate RDF for 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() print(f"Starting RDF analysis with {args.workers} workers...") main(workers=args.workers) 以上代码实现了湿法磷酸体系中水 磷酸 水合氢离子以及氟之间的RDF计算,其中将O和H分别归类。沿用该代码的框架,修改其中的判别逻辑。首先识别P,在P周围搜寻O原子,如果该O原子距离在1.6埃以内则视为Op,而对于Op在其周围搜寻H原子,如果该H在距离1.3埃以内则视为成键即该H为Hp。接着体系中全部的O原子在去除Op之后剩下的O,在这些剩余的O周围搜寻整体的H,如果H的距离在1.2埃以内则视为成键,然后依照成键的H数量判定:如果H的数量为1,则记为-OH羟基(在这里不需要计算羟基部分,只是识别出来有利于逻辑完整性,并不参与RDF计算,也不需要特别标注表明),H的数量为2,则记为H2O水(该O也随之记为Ow,对应的两个H也记为Hw),如果H的数量为3,则记为水合氢离子(该O随之记为Oh,对应的三个H也记为Hh)。体系中存在质子转移的情况,所以需要每一帧重新识别原子的归属问题,如果H同时处于两个成键识别范围则按照就近原则,离哪个近则归属到哪一个(这里包括磷酸-磷酸,磷酸-水,磷酸-水合氢离子,水-水,水-水合氢离子,水合氢离子-水合氢离子,如果H同时处于某种情况下两个化学成键范围则采用就近原则),在实时重新归属质子的情况下,计算出包含质子转移部分的RDF,在这里,我们将排除自身化学键的阈值先设置为0,不需要只看氢键部分了。
07-14
import numpy as np import matplotlib.pyplot as plt from pymatgen.io.vasp import Vasprun from pymatgen.core.structure import Structure from scipy.signal import savgol_filter from scipy.spatial import cKDTree from tqdm import tqdm import matplotlib as mpl import warnings from collections import defaultdict import os import csv import argparse import multiprocessing from functools import partial import time import dill # 忽略可能的警告 warnings.filterwarnings("ignore", category=UserWarning) # 专业绘图设置 - 符合Journal of Chemical Physics要求 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): """改进的原子类型识别函数,考虑质子转移和动态归属""" # 1. 初始化数据结构 p_oxygens = {"P=O": [], "P-O": [], "P-OH": []} phosphate_hydrogens = [] # P-OH基团中的H原子 water_oxygens = [] # 水分子中的O原子 water_hydrogens = [] # 水分子中的H原子 hydronium_oxygens = [] # 水合氢离子中的O原子 hydronium_hydrogens = [] # 水合氢离子中的H原子 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"] # 2. 构建全局KDTree用于快速查找 all_coords = np.array([site.coords for site in struct]) kdtree = cKDTree(all_coords, boxsize=struct.lattice.abc) # 3. 识别磷酸基团 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"] if not p_o_indices: continue # 按距离排序并确定P=O (最近的O) distances = [struct.get_distance(p_idx, o_idx) for o_idx in p_o_indices] sorted_indices = sorted(range(len(distances)), key=lambda i: distances[i]) closest_o = p_o_indices[sorted_indices[0]] phosphate_oxygens.append(closest_o) p_oxygens["P=O"].append(closest_o) # 处理其他O原子 for i in range(1, len(p_o_indices)): o_idx = p_o_indices[sorted_indices[i]] phosphate_oxygens.append(o_idx) p_oxygens["P-O"].append(o_idx) # 暂时标记为P-O,后续可能调整为P-OH # 4. 识别所有H原子并确定归属 hydrogen_owners = {} # 存储每个H原子的归属O原子 h_atoms = [i for i, site in enumerate(struct) if site.species_string == "H"] for h_idx in h_atoms: # 查找H周围1.2Å内的所有O原子 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 # 计算距离并找到最近的O原子 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 # 5. 重新分类磷酸氧 for o_idx in phosphate_oxygens: # 检查该氧原子是否有归属的H原子 has_hydrogen = any(owner_o == o_idx for h_idx, owner_o in hydrogen_owners.items()) # 如果是P=O且没有H,保持P=O if o_idx in p_oxygens["P=O"] and not has_hydrogen: continue # 保持为P=O # 如果是P=O但有H,需要重新分类 if o_idx in p_oxygens["P=O"] and has_hydrogen: p_oxygens["P=O"].remove(o_idx) p_oxygens["P-OH"].append(o_idx) # 如果是P-O且有H,改为P-OH if o_idx in p_oxygens["P-O"] and has_hydrogen: p_oxygens["P-O"].remove(o_idx) p_oxygens["P-OH"].append(o_idx) # 6. 识别水和水合氢离子 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原子数量 o_h_count = defaultdict(int) for h_idx, owner_o in hydrogen_owners.items(): o_h_count[owner_o] += 1 for o_idx in non_phosphate_os: h_count = o_h_count.get(o_idx, 0) # 获取归属的H原子 attached_hs = [h_idx for h_idx, owner_o in hydrogen_owners.items() if owner_o == o_idx] if h_count == 2: # 水分子 water_oxygens.append(o_idx) water_hydrogens.extend(attached_hs) elif h_count == 3: # 水合氢离子 hydronium_oxygens.append(o_idx) hydronium_hydrogens.extend(attached_hs) # h_count=1的羟基不处理 # 7. 识别磷酸基团的H原子 for o_idx in p_oxygens["P-OH"]: attached_hs = [h_idx for h_idx, owner_o in hydrogen_owners.items() if owner_o == o_idx] phosphate_hydrogens.extend(attached_hs) return { "phosphate_oxygens": p_oxygens, "phosphate_hydrogens": phosphate_hydrogens, "water_oxygens": water_oxygens, "water_hydrogens": water_hydrogens, "hydronium_oxygens": hydronium_oxygens, "hydronium_hydrogens": hydronium_hydrogens, "fluoride_atoms": fluoride_atoms, "aluminum_atoms": aluminum_atoms } # 2. RDF计算函数 - 修复负值问题和序列化问题 def process_frame(struct, center_sel, target_sel, r_max, exclude_bonds, bond_threshold): """处理单帧结构计算,完全处理空原子类型情况""" # 每帧重新识别原子类型(关键!) atom_types = identify_atom_types(struct) # 获取中心原子和目标原子 centers = center_sel(atom_types) targets = target_sel(atom_types) # 处理空原子类型情况 - 第一重保护 if len(centers) == 0 or len(targets) == 0: return { "distances": np.array([], dtype=np.float64), "n_centers": 0, "n_targets": 0, "volume": struct.volume } center_coords = np.array([struct[i].coords for i in centers]) target_coords = np.array([struct[i].coords for i in targets]) lattice = struct.lattice kdtree = cKDTree(target_coords, boxsize=lattice.abc) # 动态确定邻居数量 - 不超过目标原子数 k_val = min(50, len(targets)) # 处理目标原子数量为0的情况 - 第二重保护 if k_val == 0: return { "distances": np.array([], dtype=np.float64), "n_centers": len(centers), "n_targets": len(targets), "volume": struct.volume } # 执行查询并确保结果统一格式 try: query_result = kdtree.query(center_coords, k=k_val, distance_upper_bound=r_max) except Exception as e: # 异常处理 - 返回空结果 print(f"KDTree query error: {str(e)}") return { "distances": np.array([], dtype=np.float64), "n_centers": len(centers), "n_targets": len(targets), "volume": struct.volume } # 统一处理不同维度的返回结果 if k_val == 1: # 处理单邻居情况 if isinstance(query_result, tuple): distances, indices = query_result else: distances = query_result indices = np.zeros_like(distances, dtype=int) # 确保数组格式 distances = np.atleast_1d(distances) indices = np.atleast_1d(indices) else: # 多邻居情况 distances, indices = query_result # 确保二维数组格式 if distances.ndim == 1: distances = distances.reshape(-1, 1) indices = indices.reshape(-1, 1) valid_distances = [] for i in range(distances.shape[0]): center_idx = centers[i] for j in range(distances.shape[1]): dist = distances[i, j] # 跳过超出范围的距离 if dist > r_max or np.isinf(dist): continue target_idx = targets[indices[i, j]] # 排除化学键 if exclude_bonds: actual_dist = struct.get_distance(center_idx, target_idx) if actual_dist < bond_threshold: continue valid_distances.append(dist) return { "distances": np.array(valid_distances, dtype=np.float64), "n_centers": len(centers), "n_targets": len(targets), "volume": struct.volume } def calculate_rdf_parallel(structures, center_sel, target_sel, r_max=8.0, bin_width=0.05, exclude_bonds=True, bond_threshold=1.3, workers=1): """ 并行计算径向分布函数 :param workers: 并行工作进程数 """ bins = np.arange(0, r_max, bin_width) hist = np.zeros(len(bins) - 1) total_centers = 0 total_targets = 0 total_volume = 0 # 准备参数 - 使用dill解决序列化问题 dill.settings['recurse'] = True func = partial(process_frame, center_sel=center_sel, target_sel=target_sel, r_max=r_max, exclude_bonds=exclude_bonds, bond_threshold=bond_threshold) # 使用多进程池 with multiprocessing.Pool(processes=workers) as pool: results = [] # 使用imap_unordered提高效率 for res in tqdm(pool.imap_unordered(func, structures), total=len(structures), desc="Calculating RDF"): results.append(res) # 处理结果 - 特别注意空结果处理 n_frames = 0 for res in results: if res is None: continue n_frames += 1 valid_distances = res["distances"] n_centers = res["n_centers"] n_targets = res["n_targets"] volume = res["volume"] # 累加计数 if len(valid_distances) > 0: hist += np.histogram(valid_distances, bins=bins)[0] total_centers += n_centers total_targets += n_targets total_volume += volume # 修正归一化 - 解决负值问题 if n_frames == 0: # 没有有效帧时返回空结果 r = bins[:-1] + bin_width/2 return r, np.zeros_like(r), {"position": None, "value": None} avg_density = total_targets / total_volume if total_volume > 0 else 0 r = bins[:-1] + bin_width/2 rdf = np.zeros_like(r) for i in range(len(hist)): r_lower = bins[i] r_upper = bins[i+1] shell_vol = 4/3 * np.pi * (r_upper**3 - r_lower**3) expected_count = shell_vol * avg_density * total_centers # 避免除以零 if expected_count > 1e-10: rdf[i] = hist[i] / expected_count else: rdf[i] = 0 # 更稳健的平滑处理 - 避免边界效应 if len(rdf) > 10: window_length = min(15, len(rdf)//2*2+1) polyorder = min(5, window_length-1) rdf_smoothed = savgol_filter(rdf, window_length=window_length, polyorder=polyorder, mode='mirror') else: rdf_smoothed = rdf # 计算主要峰值 peak_info = {} mask = (r >= 1.5) & (r <= 3.0) if np.any(mask) and np.any(rdf_smoothed[mask] > 0): peak_idx = np.argmax(rdf_smoothed[mask]) peak_pos = r[mask][peak_idx] peak_val = rdf_smoothed[mask][peak_idx] peak_info = {"position": peak_pos, "value": peak_val} else: peak_info = {"position": None, "value": None} return r, rdf_smoothed, peak_info # 3. 定义精细化的选择器函数(避免lambda序列化问题) def selector_phosphate_P_double_O(atom_types): return atom_types["phosphate_oxygens"]["P=O"] def selector_phosphate_P_OH(atom_types): return atom_types["phosphate_oxygens"]["P-OH"] def selector_phosphate_P_O(atom_types): return atom_types["phosphate_oxygens"]["P-O"] def selector_phosphate_hydrogens(atom_types): return atom_types["phosphate_hydrogens"] def selector_water_only_hydrogens(atom_types): """仅选择水分子中的氢原子""" return atom_types["water_hydrogens"] def selector_hydronium_only_hydrogens(atom_types): """仅选择水合氢离子中的氢原子""" return atom_types["hydronium_hydrogens"] def selector_water_only_oxygens(atom_types): """仅选择水分子中的氧原子""" return atom_types["water_oxygens"] def selector_hydronium_only_oxygens(atom_types): """仅选择水合氢离子中的氧原子""" return atom_types["hydronium_oxygens"] def selector_fluoride_atoms(atom_types): return atom_types["fluoride_atoms"] def selector_aluminum_atoms(atom_types): return atom_types["aluminum_atoms"] def selector_all_phosphate_oxygens(atom_types): return (atom_types["phosphate_oxygens"]["P=O"] + atom_types["phosphate_oxygens"]["P-O"] + atom_types["phosphate_oxygens"]["P-OH"]) # 4. 根据您的要求定义六张图的RDF分组配置 def get_rdf_groups(): """返回六张图的RDF分组配置(完全符合您的需求)""" return { # 图1: Al的配位情况 "Al_Coordination": [ (selector_aluminum_atoms, selector_fluoride_atoms, "Al-F", "blue"), (selector_aluminum_atoms, selector_water_only_oxygens, "Al-Ow", "green"), (selector_aluminum_atoms, selector_all_phosphate_oxygens, "Al-Op", "red") ], # 图2: F与H形成的氢键 "F_Hydrogen_Bonding": [ (selector_fluoride_atoms, selector_water_only_hydrogens, "F-Hw", "lightblue"), (selector_fluoride_atoms, selector_hydronium_only_hydrogens, "F-Hh", "blue"), (selector_fluoride_atoms, selector_phosphate_hydrogens, "F-Hp", "darkblue") ], # 图3: 磷酸作为受体与周围环境的氢键(区分氧类型) "Phosphate_Acceptor": [ (selector_phosphate_P_double_O, selector_water_only_hydrogens, "P=O···Hw", "orange"), (selector_phosphate_P_double_O, selector_hydronium_only_hydrogens, "P=O···Hh", "red"), (selector_phosphate_P_O, selector_water_only_hydrogens, "P-O···Hw", "lightgreen"), (selector_phosphate_P_O, selector_hydronium_only_hydrogens, "P-O···Hh", "green"), (selector_phosphate_P_OH, selector_water_only_hydrogens, "P-OH···Hw", "lightblue"), (selector_phosphate_P_OH, selector_hydronium_only_hydrogens, "P-OH···Hh", "blue") ], # 图4: 磷酸-水-水合氢离子交叉氢键(排除同种类型) "Cross_Species_HBonding": [ (selector_phosphate_hydrogens, selector_water_only_oxygens, "Hp···Ow", "pink"), (selector_phosphate_hydrogens, selector_hydronium_only_oxygens, "Hp···Oh", "purple"), (selector_water_only_hydrogens, selector_all_phosphate_oxygens, "Hw···Op", "lightgreen"), (selector_water_only_hydrogens, selector_hydronium_only_oxygens, "Hw···Oh", "green"), (selector_hydronium_only_hydrogens, selector_water_only_oxygens, "Hh···Ow", "lightblue"), (selector_hydronium_only_hydrogens, selector_all_phosphate_oxygens, "Hh···Op", "blue") ], # 图5: 同类型分子内/间氢键(区分磷酸氧类型) "Same_Species_HBonding": [ (selector_phosphate_hydrogens, selector_phosphate_P_double_O, "Hp···P=O", "red"), (selector_phosphate_hydrogens, selector_phosphate_P_O, "Hp···P-O", "orange"), (selector_phosphate_hydrogens, selector_phosphate_P_OH, "Hp···P-OH", "yellow"), (selector_water_only_hydrogens, selector_water_only_oxygens, "Hw···Ow", "lightblue"), (selector_hydronium_only_hydrogens, selector_hydronium_only_oxygens, "Hh···Oh", "blue") ], # 图6: O-O聚集分析(Op不区分类型) "O_O_Aggregation": [ (selector_all_phosphate_oxygens, selector_water_only_oxygens, "Op-Ow", "blue"), (selector_all_phosphate_oxygens, selector_hydronium_only_oxygens, "Op-Oh", "green"), (selector_all_phosphate_oxygens, selector_all_phosphate_oxygens, "Op-Op", "red"), (selector_water_only_oxygens, selector_hydronium_only_oxygens, "Ow-Oh", "purple"), (selector_water_only_oxygens, selector_water_only_oxygens, "Ow-Ow", "cyan"), (selector_hydronium_only_oxygens, selector_hydronium_only_oxygens, "Oh-Oh", "magenta") ] } # 5. 主程序 - 优化并行处理 def main(workers=1): # 定义要处理的体系 vasprun_files = { "System1": "vasprun1.xml", "System2": "vasprun2.xml", "System3": "vasprun3.xml", "System4": "vasprun4.xml" } # 获取RDF分组配置 rdf_groups = get_rdf_groups() # 标题映射(根据您的要求) title_map = { "Al_Coordination": "Al Coordination Environment", "F_Hydrogen_Bonding": "F-H Hydrogen Bonding", "Phosphate_Acceptor": "Phosphate as H-bond Acceptor", "Cross_Species_HBonding": "Cross H-bonding between Different Species", "Same_Species_HBonding": "Intra- and Inter-molecular H-bonding", "O_O_Aggregation": "O-O Aggregation Analysis" } # 存储所有数据 all_system_data = {} group_y_max = {group_name: 0 for group_name in list(rdf_groups.keys())} group_x_max = { "Al_Coordination": (1.5, 3.5), "F_Hydrogen_Bonding": (1.0, 3.0), "Phosphate_Acceptor": (1.0, 3.0), "Cross_Species_HBonding": (1.0, 3.0), "Same_Species_HBonding": (1.0, 3.0), "O_O_Aggregation": (2.0, 6.0) } # 创建输出目录 os.makedirs("RDF_Plots", exist_ok=True) # 计算所有体系的所有RDF数据 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") # 存储体系数据 system_data = { "rdf_results": {}, "peak_infos": {} } # 计算所有RDF分组 for group_name, pairs in rdf_groups.items(): system_data["rdf_results"][group_name] = {} system_data["peak_infos"][group_name] = {} group_y_max_current = 0 for center_sel, target_sel, label, color in pairs: print(f"\nCalculating RDF for: {label}") try: r, rdf, peak_info = calculate_rdf_parallel( structures, center_sel, target_sel, r_max=10.0, exclude_bonds=True, bond_threshold=1.3, workers=workers ) system_data["rdf_results"][group_name][label] = (r, rdf, color) system_data["peak_infos"][group_name][label] = peak_info if len(rdf) > 0: current_max = np.max(rdf) if current_max > group_y_max_current: group_y_max_current = current_max if peak_info["position"] is not None: print(f" Peak for {label}: {peak_info['position']:.3f} Å (g(r) = {peak_info['value']:.2f})") else: print(f" No significant peak found for {label} in 1.5-3.0 Å range") except Exception as e: print(f"Error calculating RDF for {label}: {str(e)}") system_data["rdf_results"][group_name][label] = (np.array([]), np.array([]), color) system_data["peak_infos"][group_name][label] = {"position": None, "value": None} if group_y_max_current > group_y_max[group_name]: group_y_max[group_name] = group_y_max_current all_system_data[system_name] = system_data 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)}") # 为每个分组添加余量 for group_name in group_y_max: group_y_max[group_name] = max(group_y_max[group_name] * 1.15, 3.0) # 确保最小值 # 第二步:生成符合期刊要求的图表 for system_name, system_data in all_system_data.items(): print(f"\nGenerating publication-quality plots for {system_name}") for group_name, group_data in system_data["rdf_results"].items(): fig, ax = plt.subplots(figsize=(8, 6)) # 设置坐标轴范围 xlim = group_x_max.get(group_name, (0, 6.0)) ylim = (0, group_y_max[group_name]) for label, (r, rdf, color) in group_data.items(): if len(r) > 0 and len(rdf) > 0: ax.plot(r, rdf, color=color, label=label, linewidth=2.0) ax.set_xlim(xlim) ax.set_ylim(ylim) # 期刊格式标签 ax.set_xlabel('Radial Distance (Å)', fontweight='bold') ax.set_ylabel('g(r)', fontweight='bold') # 添加体系名称到标题 ax.set_title(f"{system_name}: {title_map[group_name]}", fontsize=16, pad=15) # 精简图例 ncol = 3 if group_name == "Same_Species_HBonding" else 1 # 图5使用三列图例 ax.legend(ncol=ncol, loc='best', framealpha=0.8, fontsize=10) # 添加氢键区域标记(除O-O聚集图外) if group_name != "O_O_Aggregation": ax.axvspan(1.5, 2.5, alpha=0.1, color='green', zorder=0) ax.text(1.7, ylim[1]*0.85, 'H-bond Region', fontsize=10) # 添加网格 ax.grid(True, linestyle='--', alpha=0.5) # 保存高分辨率图片 plt.tight_layout() filename = os.path.join("RDF_Plots", f"RDF_{system_name}_{group_name}.tiff") plt.savefig(filename, bbox_inches='tight', dpi=600, format='tiff') print(f"Saved publication plot: {filename}") plt.close() # 保存Origin兼容数据 save_origin_data(system_name, system_data) print("\nAll RDF analysis completed successfully!") def save_origin_data(system_name, system_data): """保存Origin兼容格式数据""" os.makedirs("Origin_Data", exist_ok=True) system_dir = os.path.join("Origin_Data", system_name) os.makedirs(system_dir, exist_ok=True) # 保存峰值信息 peak_info_path = os.path.join(system_dir, f"Peak_Positions_{system_name}.csv") with open(peak_info_path, 'w', newline='') as csvfile: writer = csv.writer(csvfile) writer.writerow(["Group", "Interaction", "Peak Position (A)", "g(r) Value"]) for group_name, peaks in system_data["peak_infos"].items(): for label, info in peaks.items(): if info["position"] is not None: writer.writerow([group_name, label, f"{info['position']:.3f}", f"{info['value']:.3f}"]) else: writer.writerow([group_name, label, "N/A", "N/A"]) print(f"Saved peak positions: {peak_info_path}") # 保存RDF数据 for group_name, group_results in system_data["rdf_results"].items(): group_dir = os.path.join(system_dir, group_name) os.makedirs(group_dir, exist_ok=True) for label, (r, rdf, color) in group_results.items(): if len(r) > 0 and len(rdf) > 0: safe_label = label.replace(" ", "_").replace("/", "_").replace("=", "_") safe_label = safe_label.replace("(", "").replace(")", "").replace("$", "") filename = f"RDF_{system_name}_{group_name}_{safe_label}.csv" filepath = os.path.join(group_dir, filename) with open(filepath, 'w', newline='') as csvfile: writer = csv.writer(csvfile) writer.writerow(["Distance (A)", "g(r)"]) for i in range(len(r)): writer.writerow([f"{r[i]:.6f}", f"{rdf[i]:.6f}"]) print(f"Saved Origin data: {filename}") if __name__ == "__main__": # 设置命令行参数 parser = argparse.ArgumentParser(description='Calculate RDF for 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() print(f"Starting RDF analysis with {args.workers} workers...") main(workers=args.workers) 以上代码计算的VASP数据,同时进行了详细区分计算,以及计算条件,在这里我们利用同样的框架,但修改其内容,在这里我们只需要计算体系中存在的各个元素(可能换体系,元素可能会涉及Al Fe Mg F H O P Si)之间的RDF计算,所以在这里代码当中可能需要添加可更改的中心原子和目标原子,并输出文本结果
07-15
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值