LAMMPS完成reaxff模拟后统计产物分子

输入脚本示例

# .....
units           real
atom_style      charge
read_data       data.CHO
pair_style      reax/c NULL
pair_coeff      * * ffield.reax.cho H C O 
neighbor        2 bin 
neigh_modify    every 10 delay 0 check no
fix             1 all nvt temp 3500.0 3500.0 100.0
fix             2 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
fix             4 all reax/c/species 1 100 100 species.out element H C O cutoff 2 2 0.15 
velocity        all create 3500 428459 dist gaussian
timestep        0.25
thermo          1  
run             3000

注：velocity命令给体系赋予一个初始速度，不加这条命令则体系从0K开始升温到3500K。428459是种子号，改成12345也可以。

统计产物的核心命令

fix reax/species
fix 1 all reax/species 1 100 100 species.out element H C O cutoff 2 2 0.15

1 100 100 分别表示每隔1步采样一次键级，采样100次，并在100步的时候输出平均键级，再比如改成1 2 100则表示在98 100步时采集键级在100步时取平均输出。species.out则是输出文件的文件名。element参数建议加上，不过不加则默认为C H O的顺序，而data文件中的原子顺序大多数情况下并非如此。cutoff也是可选参数，表示原子之间的截止键级（注意是键级不是距离），需要自己酌情调整，数值太大会导致完整的分子被切成碎片，数值太小则会导致断裂的分子之间藕断丝连。

运行LAMMPS后species.out文件格式

# Timestep     No_Moles     No_Specs     H11C5  H  H8C5  H2  H7C5  H4C  H4C3  H3C  H9C5  O2  100         29         10   1   10   1   2   1   1   1   1   1   10  
# Timestep     No_Moles     No_Specs     H5C2  H5C3  H2  H8C5  H  H4C3  H3C2  H3C  H2C  H4C2  O2  200         33         11   1   1   5   1   7   2   2   1   2   1   10  # Timestep     No_Moles     No_Specs     H5C2  H3C2  H2C  H2  H4C3  H  H4C2  H3C  HC  O2  300         36         10   1   3   4   6   3   6   1   1   1   10  
# Timestep     No_Moles     No_Specs     H5C2  H3C2  H2C  H2  H4C3  H  H4C2  H3C  HC  O2  400         36         10   1   3   4   6   3   6   1   1   1   10  
# Timestep     No_Moles     No_Specs     H5C2  H3C2  H  H2C  H2  H4C3  H4C2  H3C  HC  O2  500         36         10   1   3   6   4   6   3   1   1   1   10  ... ...

.out文件中奇数行为注释行，标注了一些分子、团簇的符号，偶数行则是团簇分子的个数。

一个简单的python脚本读取指定分子的个数。

# Setting parameters
file_name='pecies.out'
target='O2'
# Read coordinate from dump file 
z=[]
type=[]
with open(file_name,'r') as species:
    lines=species.readlines()
    length=len(lines) 
species.close()
Step=[]
Num=[]
for i in range(0,length,2):
    a=lines[i].strip()    
    a=a.split()    
    b=lines[i+1].strip()    
    b=b.split() 
    if target in a:
        Num.append(b[a.index(target)-1]) 
    else:
        Num.append(0)
        Step.append(b[0])
print('#Time_step  ',target)
for i in range(len(Num)): 
    print (Step[i],Num[i])