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计算,所以在这里代码当中可能需要添加可更改的中心原子和目标原子,并输出文本结果