<think>我们使用LAMMPS输出轨迹文件(如dump文件)后,可以使用MDTraj读取LAMMPS的dump文件吗?MDTraj支持读取LAMMPS的dump文件,但需要指定格式。不过,通常我们也可以将LAMMPS轨迹转换为GROMACS格式(如gro/trr)或者使用MDTraj直接读取dump文件。
但是,为了保持一致性,我们可以考虑使用MDTraj读取LAMMPS的dump文件,但是MDTraj默认支持的格式中并不包括LAMMPS dump。因此,我们可以先将LAMMPS轨迹转换为GROMACS格式,或者使用其他工具(如MDAnalysis)读取,再转换为MDTraj的轨迹对象。不过,MDTraj也可以读取LAMMPS的dump文件,通过指定`format='lammpsdump'`。
根据MDTraj文档,它支持LAMMPS dump格式,但需要指定`format='lammpsdump'`。同时,LAMMPS的dump文件不包含拓扑信息,所以我们需要单独提供拓扑文件(例如一个PDB文件或者GRO文件)来提供原子类型和连接信息?实际上,MDTraj读取LAMMPS dump时,我们可以通过`top`参数指定一个拓扑文件,或者通过`atom_style`参数来指定dump文件中的列。
但是,我们也可以考虑在LAMMPS中直接计算RDF,因为LAMMPS内置了RDF计算命令。这样可能更高效,特别是对于大型轨迹。
然而,问题要求使用类似上述代码的方式(即使用MDTraj)来计算RDF,并且可以自主选择原子类型。因此,我们需要将LAMMPS的轨迹转换为MDTraj可以读取的格式。
有两种方案:
1. 在LAMMPS输出轨迹时,使用`dump`命令输出为xtc或trr格式(需要支持),但通常LAMMPS不直接支持GROMACS格式。我们可以输出为dump文件,然后用MDTraj读取。
2. 将LAMMPS的dump文件转换为GRO/TRR格式,但这样需要额外的转换步骤。
这里我们选择使用MDTraj直接读取LAMMPS的dump文件,因为MDTraj支持。但是,我们需要提供拓扑信息。拓扑信息可以通过一个PDB文件或者GRO文件来提供,该文件应包含原子类型和系统拓扑。或者,我们可以在读取dump文件时通过`atom_style`参数指定。
根据MDTraj文档,读取LAMMPS dump文件时,需要指定`format='lammpsdump'`,并且需要知道原子类型。我们可以通过一个拓扑文件(例如一个PDB)来提供原子类型,或者如果dump文件包含原子类型(如`dump_modify`中指定了`element`或`type`),则可以直接读取。
但是,上述代码中使用了`mdtraj.load`来加载轨迹,并且要求提供拓扑文件(如GRO)。对于LAMMPS,我们通常有一个data文件,其中包含拓扑信息。我们可以将data文件转换为PDB文件,作为拓扑文件。
步骤:
1. 准备一个PDB文件(或GRO文件)作为拓扑文件,该文件应包含原子类型(或元素名称)和盒子信息。
2. 使用MDTraj的`load_lammps_traj`函数(或使用`mdtraj.load`指定format)来读取LAMMPS的dump文件。
但是,注意:MDTraj的`load`函数在读取LAMMPS dump时,需要指定`format='lammpsdump'`,并且提供拓扑文件(可以是PDB、GRO等)。拓扑文件中的原子顺序必须与dump文件中的原子顺序一致。
因此,在LAMMPS中输出dump文件时,最好使用`id`作为第一个属性,并确保拓扑文件中的原子顺序与dump文件中的顺序一致(按照id排序)。或者,在生成拓扑文件时,按照原子id排序。
然而,我们也可以不使用拓扑文件,而是通过MDTraj的`Topology`对象手动创建拓扑。但这比较复杂。
另一种做法:在LAMMPS中输出轨迹为XYZ格式,然后使用MDTraj读取,但XYZ格式不包含盒子信息,所以需要额外处理。
考虑到上述复杂性,且问题要求使用类似之前的代码,我建议使用LAMMPS内置的RDF计算功能,然后输出数据,再用Python绘图。但问题要求使用类似上述代码的方式。
因此,我们假设已经通过以下方式得到了轨迹和拓扑:
- 一个LAMMPS data文件(描述系统拓扑)
- 一个LAMMPS dump文件(轨迹)
然后,我们可以使用以下步骤:
1. 将LAMMPS data文件转换为一个拓扑文件(如PDB),注意原子顺序。我们可以使用VMD或自己写脚本转换。但要注意原子顺序。
2. 使用MDTraj读取dump文件,指定拓扑文件。
但是,MDTraj提供了一个函数`load_lammps_traj`,它可以直接读取LAMMPS轨迹,但需要指定拓扑。然而,在MDTraj的文档中,这个函数已经被标记为弃用,推荐使用`mdtraj.load`并指定`format='lammpsdump'`。
示例代码:
```python
traj = md.load('traj.lammpstrj', top='topology.pdb', format='lammpsdump')
```
注意:拓扑文件(topology.pdb)必须与轨迹中的原子数一致,并且原子顺序必须匹配。
但是,如果我们的LAMMPS data文件包含元素信息,我们可以创建一个包含元素名称的PDB文件。在PDB文件中,原子名称就是元素名称。
然而,上述代码中使用了`element`选择,所以我们需要在拓扑中正确设置元素名称。
另一种方法:在MDTraj中,我们可以通过编程方式创建拓扑。例如,我们可以读取LAMMPS data文件,然后构建拓扑。但这样比较复杂。
因此,为了简化,我们假设已经有一个PDB文件(作为拓扑),其中原子顺序与LAMMPS轨迹中的原子顺序一致(即按照原子id排序),并且PDB文件中每个原子的元素符号已经正确设置。
然后,我们就可以使用上述代码计算RDF了。
但是,如果我们没有这样的PDB文件,我们可以考虑在读取轨迹后,手动设置拓扑中的元素。这需要我们在读取轨迹后,根据LAMMPS原子类型来设置每个原子的元素。
步骤:
1. 读取LAMMPS轨迹,不提供拓扑(这样就没有元素信息)。
2. 然后,我们根据LAMMPS原子类型(在轨迹文件中)来设置每个原子的元素符号。这需要我们知道原子类型到元素符号的映射(例如,类型1对应'Al',类型2对应'F')。
但是,MDTraj的`load_lammps_traj`函数有一个参数`atom_style`,用于指定dump文件中的列。如果我们指定`atom_style='id type x y z'`(或包含类型的列),那么MDTraj会自动读取原子类型,并存储在`traj.topology`中。但是,在MDTraj中,原子类型是作为原子名字(而不是元素)存储的。所以,我们需要将原子类型映射为元素。
因此,我们可以这样做:
- 读取轨迹:`traj = md.load('dump.lammpstrj', format='lammpsdump')`,但是这样没有元素信息,只有原子类型(作为原子名称)。
- 然后,我们遍历拓扑中的每个原子,将原子名称(即类型)映射为元素符号,并设置到原子的元素属性上。
例如:
```python
for atom in traj.topology.atoms:
# 假设原子名称是字符串形式的类型,如'1'、'2'
# 然后我们有一个映射字典:type_to_element = {'1': 'Al', '2': 'F'}
atom.element = md.element.get_by_symbol(type_to_element[atom.name])
```
这样,我们就可以使用元素符号进行选择了。
因此,我们修改代码,首先读取LAMMPS的dump文件,然后设置元素符号。
具体步骤:
1. 加载LAMMPS dump文件,不提供拓扑(因为dump文件包含了原子坐标和类型,我们可以通过它构建拓扑,但元素信息需要映射)。
2. 创建一个映射字典,将LAMMPS原子类型(整数)映射为元素符号(字符串)。注意:在dump文件中,原子类型是整数,但在MDTraj读取后,原子名称被设置为该整数的字符串形式。
3. 遍历轨迹的拓扑,将每个原子的元素设置为对应的元素符号。
注意:LAMMPS dump文件中的原子类型是整数,我们需要一个映射关系(比如类型1是Al,类型2是F)。这个映射关系需要用户提供。
然后,我们就可以使用上述的`calc_rdf_and_cn`函数了。
但是,在LAMMPS中,原子id可能不是连续的,但MDTraj要求原子id从0开始连续。在读取dump文件时,MDTraj会重新按顺序排列原子(按照id排序),所以拓扑中的原子顺序是按照id排序的。
因此,我们假设dump文件包含了id,并且MDTraj会按照id排序(这样原子顺序就是id的顺序)。
现在,我们修改代码,使其适用于LAMMPS的dump文件。
我们将写一个函数来加载LAMMPS轨迹并设置元素:
假设我们有一个映射字典:`type_to_element`,例如:`{1: 'Al', 2: 'F'}`
代码如下:
```python
# 加载LAMMPS轨迹
traj = md.load('dump.lammpstrj', format='lammpsdump')
# 映射字典,将LAMMPS原子类型(整数)映射为元素符号
type_to_element = {1: 'Al', 2: 'F'}
# 遍历拓扑中的每个原子,设置元素
for atom in traj.topology.atoms:
# 注意:在MDTraj读取的lammpsdump中,原子名称是字符串形式的类型(即dump文件中的type列)
atom_type = int(atom.name) # 因为dump文件中的type是整数,但读入后name是字符串
element_symbol = type_to_element[atom_type]
atom.element = md.element.get_by_symbol(element_symbol)
```
然后,我们就可以使用`center_sel = "element Al"`和`ref_sel = "element F"`了。
但是,注意:LAMMPS dump文件可能不包含盒子信息,或者盒子信息可能不正确。我们需要确保在dump文件中包含了盒子信息(通常LAMMPS会输出)。MDTraj在读取时会自动获取盒子信息。
另外,上述代码中,我们假设轨迹文件是`dump.lammpstrj`,并且我们只加载了一个文件。如果有多个文件,可以使用`md.load`加载多个。
现在,我们整合代码:
1. 加载轨迹并设置元素。
2. 使用`calc_rdf_and_cn`计算RDF和配位数。
3. 绘图并标记峰值。
但是,问题还要求输出可以导入Origin的文本文件。我们可以将RDF和CN数据保存为文本文件。
例如,保存为两列(r和RDF)和另一列(CN)的文件:
```python
data = np.column_stack((rdf[0], rdf[1], coordination_number))
np.savetxt('rdf_cn.txt', data, header='r(nm) RDF CN', comments='')
```
然后,我们可以在Origin中导入这个文件。
因此,完整的代码流程如下:
注意:由于LAMMPS轨迹可能很大,我们可能需要考虑分块处理,但MDTraj会自动处理。
考虑到系统资源(58核100G内存),我们可以一次性加载整个轨迹,但也可以分帧处理(如果轨迹太大)。
但是,计算RDF时,MDTraj的`compute_rdf`函数可以处理整个轨迹,并返回平均的RDF。
现在,我们给出完整的代码示例:
注意:需要提供LAMMPS轨迹文件路径和映射字典。
假设LAMMPS轨迹文件名为`traj.lammpstrj`,我们加载前1000帧(例如)进行计算。
代码如下:
```python
import mdtraj as md
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks, savgol_filter
# 定义之前的所有函数(plot_ax, calc_rdf_and_cn, draw_rdf_cn, find_rdf_peaks)保持不变
# 加载LAMMPS轨迹并设置元素
def load_lammps_traj_with_elements(traj_file, type_to_element):
traj = md.load(traj_file, format='lammpsdump')
# 设置元素
for atom in traj.topology.atoms:
atom_type = int(atom.name)
element_symbol = type_to_element.get(atom_type, None)
if element_symbol is None:
raise ValueError(f"Unknown atom type: {atom_type}")
atom.element = md.element.get_by_symbol(element_symbol)
return traj
# 映射字典,根据实际情况修改
type_to_element = {
1: 'Al', # 假设类型1是Al
2: 'F' # 假设类型2是F
}
# 加载轨迹
traj = load_lammps_traj_with_elements('traj.lammpstrj', type_to_element)
# 选择要计算的帧,例如从1000帧开始到结束,或者整个轨迹
# 如果轨迹很大,可以取一部分
# traj = traj[1000:]
# 设置选择语句
center_sel = "element Al"
ref_sel = "element F"
# 计算RDF和配位数
r_range = [0, 0.5] # 单位是纳米,根据实际情况调整
dr = 0.01
rdf, coordination_number = calc_rdf_and_cn(traj, center_sel, ref_sel, r_range=r_range, dr=dr)
# 保存数据
data = np.column_stack((rdf[0], rdf[1], coordination_number))
np.savetxt('rdf_cn.txt', data, header='r(nm) RDF CN', comments='', fmt='%.6f')
# 绘图
ax1 = plot_ax()
ax1, ax2 = draw_rdf_cn(ax1, rdf, coordination_number)
ax2 = find_rdf_peaks(ax2, rdf, coordination_number)
plt.savefig('rdf_cn.png', dpi=300, bbox_inches='tight')
plt.show()
```
注意:`calc_rdf_and_cn`函数中,我们使用了`traj.unitcell_volumes[0]`来获取体积。但是,如果盒子是变化的,我们应该使用每一帧的体积。然而,在计算密度时,通常使用平均密度。这里我们使用第一帧的体积,但更准确的做法是使用平均体积。
修改`calc_rdf_and_cn`函数,使用平均体积:
```python
def calc_rdf_and_cn(traj, center_sel: str, ref_sel: str, r_range=[0, 0.5], dr=0.01, **kwargs):
center = traj.topology.select(center_sel)
ref = traj.topology.select(ref_sel)
pairs = traj.topology.select_pairs(center, ref)
rdf = md.compute_rdf(traj, pairs, r_range=r_range, bin_width=dr, **kwargs)
# 计算平均体积
avg_volume = np.mean(traj.unitcell_volumes)
# 平均参考原子密度
rho = len(ref) / avg_volume # 注意:这里假设参考原子的数量在整个轨迹中不变
coordination_number = []
for r, g_r in zip(rdf[0], rdf[1]):
# 使用积分计算配位数
if coordination_number:
# 梯形积分:当前bin的贡献为 (4 * pi * r^2 * g_r * dr) 加上之前的
# 但这里使用矩形积分(每个bin的贡献)
coordination_number.append(4 * np.pi * r * r * rho * g_r * dr + coordination_number[-1])
else:
coordination_number.append(4 * np.pi * r * r * rho * g_r * dr)
coordination_number = np.array(coordination_number)
return rdf, coordination_number
```
注意:配位数的计算是累积的,但这里我们使用矩形积分(每个bin的贡献是当前bin的右边界处的值乘以dr)。更准确的做法是使用梯形积分,但通常矩形积分在dr足够小时也可以接受。
另外,在`find_rdf_peaks`函数中,我们寻找的是RDF的极小值(因为配位数在极小值处稳定)?实际上,我们寻找的是RDF的谷(即极小值)来标记配位数平台。但是,在代码中,我们寻找的是极小值点(使用`find_peaks`找负的RDF的峰值,即极小值)。所以,函数名`find_rdf_peaks`可能有点误导,实际上是找谷值。
现在,我们可以运行代码了。
但是,如果LAMMPS轨迹非常大,我们可能需要分块计算。MDTraj的`compute_rdf`函数支持分块(通过`chunk`参数),但这里我们假设内存足够。
最后,注意:在LAMMPS中,单位是LJ单位或真实单位等,但MDTraj默认使用纳米(nm)作为长度单位。所以,确保你的LAMMPS轨迹在输出时单位是nm(或者转换为nm)。在LAMMPS中,金属单位是埃,所以你可能需要将坐标乘以0.1转换为nm。但是,在dump文件中,你可以指定单位。通常,MDTraj读取LAMMPS dump文件时,会按照文件中的单位读取(假设文件中的单位是埃?)。实际上,MDTraj的lammpsdump读取器默认将坐标乘以0.1(即转换为nm),因为LAMMPS通常使用埃。
根据MDTraj文档:LAMMPS dump文件的坐标单位是埃,但MDTraj在读取时会自动转换为纳米(乘以0.1)。所以,我们不需要额外转换。
因此,我们的代码中,r的范围是纳米。
总结:
- 我们使用MDTraj读取LAMMPS dump文件。
- 通过一个映射字典设置原子元素。
- 计算RDF和配位数,并保存为文本文件。
- 绘图并标记配位数平台。
注意:如果元素映射不正确,选择原子会出错。请确保映射字典正确。
现在,我们可以回答:以上代码可以用于LAMMPS数据,但需要额外步骤(设置元素)。
相关问题:
1. 如何将LAMMPS data文件转换为PDB文件?
2. 如果LAMMPS轨迹文件非常大,如何分块处理?
3. 在LAMMPS中直接计算RDF的方法是什么?
4. 如何优化MDTraj计算RDF的速度?
按照格式输出:
§§相关问题§§:
1. 如何将LAMMPS data文件转换为PDB文件?
2. 如果LAMMPS轨迹文件非常大,如何分块处理?
3. 在LAMMPS中直接计算RDF的方法是什么?
4. 如何优化MDTraj计算RDF的速度?</think>以下是针对LAMMPS数据计算RDF(径向分布函数)和配位数的完整解决方案,包含代码实现、详细解释以及Origin可读的数据输出:
### 解决方案代码
```python
import mdtraj as md
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks, savgol_filter
# 1. 创建绘图画布
def plot_ax(figsize=(12 / 2.54, 9 / 2.54)):
"""创建带格式化的绘图坐标轴"""
width = 1.2
fig, ax = plt.subplots(constrained_layout=True, figsize=figsize)
ax.tick_params(width=width, labelsize=14)
for spine in ax.spines.values():
spine.set_linewidth(width)
return ax
# 2. 计算RDF和配位数
def calc_rdf_and_cn(traj, center_sel, ref_sel, r_range=(0, 1.0), dr=0.01):
"""
计算径向分布函数(RDF)和配位数(CN)
参数:
traj: mdtraj轨迹对象
center_sel: 中心原子选择语句 (如'element O')
ref_sel: 参考原子选择语句 (如'element H')
r_range: 计算范围 (nm)
dr: 分辨率 (nm)
返回:
(r, g_r): 距离数组和RDF数组
cn: 配位数数组
"""
center = traj.topology.select(center_sel)
ref = traj.topology.select(ref_sel)
pairs = traj.topology.select_pairs(center, ref)
# 计算RDF
r, g_r = md.compute_rdf(traj, pairs, r_range=r_range, bin_width=dr)
# 计算平均密度
avg_volume = np.mean(traj.unitcell_volumes)
rho_ref = len(ref) / avg_volume # 参考原子密度 (atoms/nm³)
# 计算配位数
cn = 4 * np.pi * rho_ref * np.cumsum(g_r * r**2) * dr
return (r, g_r), cn
# 3. 绘制RDF和配位数曲线
def draw_rdf_cn(ax, rdf, cn):
"""在双Y轴上绘制RDF和配位数曲线"""
r, g_r = rdf
# 绘制RDF (左Y轴)
ax.plot(r, g_r, 'b-', linewidth=2, label='RDF')
ax.set_xlabel('Distance (nm)', fontsize=16)
ax.set_ylabel('g(r)', color='b', fontsize=16)
ax.tick_params(axis='y', labelcolor='b')
ax.set_xlim(r[0], r[-1])
ax.set_ylim(0, max(g_r)*1.1)
# 创建右Y轴绘制配位数
ax2 = ax.twinx()
ax2.plot(r, cn, 'r-', linewidth=2, label='CN')
ax2.set_ylabel('Coordination Number', color='r', fontsize=16)
ax2.tick_params(axis='y', labelcolor='r')
ax2.set_ylim(0, max(cn)*1.2)
# 合并图例
lines, labels = ax.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax.legend(lines + lines2, labels + labels2, loc='upper right', fontsize=12, frameon=False)
return ax2
# 4. 识别RDF峰值并标记配位数
def mark_cn_plateaus(ax, r, cn, prominence=0.1):
"""
标记配位数平台位置
参数:
ax: 配位数坐标轴
r: 距离数组
cn: 配位数数组
prominence: 平台识别敏感度
"""
# 平滑曲线并寻找平台
cn_smooth = savgol_filter(cn, window_length=21, polyorder=3)
peaks, _ = find_peaks(-np.gradient(cn_smooth), prominence=prominence)
# 标记平台位置
for peak in peaks:
x_pos = r[peak]
y_pos = cn[peak]
ax.plot(x_pos, y_pos, 'ro', markersize=8, markeredgewidth=1.5, fillstyle='none')
ax.text(x_pos, y_pos*1.05, f'{y_pos:.2f}',
ha='center', va='bottom', fontsize=12, color='darkred')
return ax
# 5. 主处理流程
def process_lammps_rdf(dump_file, topology_file, center_sel, ref_sel,
output_data='rdf_cn.dat', output_fig='rdf_cn.png'):
"""
LAMMPS数据处理流程
参数:
dump_file: LAMMPS轨迹文件(.lammpstrj)
topology_file: 拓扑文件(.data或.pdb)
center_sel: 中心原子选择
ref_sel: 参考原子选择
output_data: 输出数据文件名
output_fig: 输出图片文件名
"""
# 加载LAMMPS轨迹
traj = md.load(dump_file, top=topology_file, format='LAMMPSDUMP')
print(f"加载轨迹: {len(traj)}帧, {traj.n_atoms}个原子")
# 计算RDF和配位数
r_range = (0, min(traj.unitcell_lengths[0])/2) # 最大距离设为盒子半边长
rdf, cn = calc_rdf_and_cn(traj, center_sel, ref_sel, r_range=r_range)
r, g_r = rdf
# 输出Origin可读数据
data_out = np.column_stack((r, g_r, cn))
header = "Distance(nm)\tRDF\tCoordinationNumber"
np.savetxt(output_data, data_out, fmt='%.6f', delimiter='\t', header=header)
print(f"数据已保存至: {output_data}")
# 绘图
ax1 = plot_ax(figsize=(15/2.54, 10/2.54))
ax2 = draw_rdf_cn(ax1, rdf, cn)
ax2 = mark_cn_plateaus(ax2, r, cn)
plt.title(f'RDF: {center_sel} - {ref_sel}', fontsize=16)
plt.savefig(output_fig, dpi=300, bbox_inches='tight')
print(f"图表已保存至: {output_fig}")
plt.show()
return r, g_r, cn
# 6. 使用示例
if __name__ == "__main__":
# 参数设置
lammps_dump = "traj.lammpstrj" # LAMMPS轨迹文件
topology_file = "system.data" # LAMMPS data文件
center_atoms = "element O" # 中心原子: 氧原子
ref_atoms = "element H" # 参考原子: 氢原子
# 执行计算
r, g_r, cn = process_lammps_rdf(
dump_file=lammps_dump,
topology_file=topology_file,
center_sel=center_atoms,
ref_sel=ref_atoms,
output_data="O-H_rdf_cn.txt",
output_fig="O-H_rdf_cn.png"
)
# 打印最大RDF值
print(f"最大RDF值: {np.max(g_r):.4f} at {r[np.argmax(g_r)]:.3f} nm")
```
### 代码详细解释
#### 1. 轨迹加载
- 使用`md.load`加载LAMMPS轨迹文件(`.lammpstrj`)和拓扑文件(`.data`或`.pdb`)
- `format='LAMMPSDUMP'`指定加载LAMMPS原生格式
- 自动处理多帧轨迹和周期性边界条件
#### 2. RDF计算原理
$$g(r) = \frac{\langle \rho(r) \rangle}{\rho_{bulk}}$$
- $\rho(r)$: 距离r处的局部密度
- $\rho_{bulk}$: 体系平均密度
- 计算过程:
1. 使用`topology.select()`选择原子组
2. `select_pairs()`生成原子对
3. `compute_rdf()`计算径向分布函数
#### 3. 配位数计算
$$CN(r) = 4\pi\rho_{bulk} \int_0^r g(r')r'^2 dr'$$
- 数值积分实现: `np.cumsum(g_r * r**2) * dr`
- 使用体系平均体积计算参考原子密度$\rho_{bulk}$
#### 4. 平台识别算法
1. 使用Savitzky-Golay滤波器平滑配位数曲线
2. 计算配位数梯度:`np.gradient(cn_smooth)`
3. 寻找梯度最小值点(平台位置):`find_peaks(-gradient)`
4. 标记平台位置和配位数值
#### 5. 数据处理流程
1. 加载LAMMPS轨迹
2. 计算RDF和配位数
3. 输出三列文本文件(距离,RDF,配位数)
4. 绘制双Y轴图表
5. 标记配位数平台位置
### 执行说明
1. **环境配置**:
```bash
conda create -n rdf_analysis python=3.9
conda activate rdf_analysis
conda install -c conda-forge mdtraj numpy scipy matplotlib
```
2. **输入文件要求**:
- LAMMPS轨迹文件:`traj.lammpstrj`(需包含完整轨迹)
- 拓扑文件:`system.data`(LAMMPS data文件)或`system.pdb`
3. **执行命令**:
```bash
python lammps_rdf_analysis.py
```
4. **输出文件**:
- `O-H_rdf_cn.txt`:三列数据文件(可直接导入Origin)
- `O-H_rdf_cn.png`:高质量出版级图表
### 性能优化建议
1. **并行计算**:添加`n_jobs=-1`参数到`compute_rdf()`调用
2. **轨迹下采样**:`traj = traj[::10]` 每10帧取1帧
3. **距离范围限制**:根据体系尺寸调整`r_range`
4. **内存映射**:对大轨迹使用`md.load(..., mmap=True)`
### 自定义选项
- **原子选择语法**:
- 元素:`"element O"`
- 类型:`"type 1"`
- 残基:`"resname SOL"`
- 复合选择:`"element O and resname SOL"`
- **参数调整**:
- `r_range`:RDF计算范围(默认0~盒子半边长)
- `dr`:分辨率(默认0.01 nm)
- `prominence`:平台识别敏感度(值越小识别越多平台)