in文件——1.ave/chunk计算层间分布情况、密度分布函数 2. 径向分布函数RDF

1.ave/chunk计算层间分布情况

#log mmtnpt.log
variable filename1 index density

# ------------------------ INITIALIZATION ----------------------------
echo screen
dimension 3
boundary p p p

atom_style full
units real

# ----------------------- ATOM DEFINITION ----------------------------
pair_style lj/cut/coul/long 12.5 8.5
pair_modify mix arithmetic #forcefield see in data file,系数组合方式
bond_style harmonic
kspace_style ewald 1.0e-6

read_restart equilnpt100ps.restart

neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes

# ----------------------- Equilibration ----------------------------
variable T0 equal 298.0 #初始温度设定为298.0K（常温）
variable T1 equal 298.0 #终止温度设定为298.0K（常温）
variable P0 equal 1.0 #初始压强设定为1个atm（常压），0.1MPa=100KPa
variable P1 equal 1.0 #终止压强设定为1个atm（常压）
variable n equal count(all) #计算所有原子的个数
variable p equal etotal/$n

################################# 1. group ############################
group cation type 5 # Na 原子，层间阳离子
group oxygen type 12 # O 原子
group hydrogen type 13 # H 原子
group water type 12 13 # water

################################# 3. 系综计算 ############################
dump 3 all custom 10 ${filename1}.lammpstrj id type x y z vx vy vz fx fy fz #可看看平衡的轨迹

reset_timestep 0
timestep 1
fix 3 all npt temp ${T0} ${T1} 100 iso ${P0} ${P1} 1000 drag 1 #温度系数步长100倍，压强系数取步长1000倍 ,各向同性控压0 atm
thermo 10
thermo_style custom step temp vol lx ly lz press pxx pyy pzz pe ke etotal

################################# 2. 分chunk计算 ############################
compute 1 cation chunk/atom bin/1d z 0 0.05
fix output_cation cation ave/chunk 10 10 100 1 density/mass density/number vx vy vz file chunk_cation.dat

compute 2 oxygen chunk/atom bin/1d z 0 0.05
fix output_oxygen oxygen ave/chunk 10 10 100 2 density/mass density/number vx vy vz file chunk_oxygen.dat

compute 3 hydrogen chunk/atom bin/1d z 0 0.05
fix output_hydrogen hydrogen ave/chunk 10 10 100 3 density/mass density/number vx vy vz file chunk_hydrogen.dat

run 100

write_restart ${filename1}.restart
write_data ${filename1}.data

 

 

2.RDF

variable filename1 index 01rdf-new

# ------------------------ INITIALIZATION ----------------------------
echo screen
dimension 3
boundary p p p

atom_style full
units real

# ----------------------- ATOM DEFINITION ----------------------------
pair_style lj/cut/coul/long 12.5 8.5
pair_modify mix arithmetic #forcefield see in data file,系数组合方式
bond_style harmonic
kspace_style ewald 1.0e-6

read_restart equilnpt100ps.restart #读入npt平衡后导出来的文件

neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes

reset_timestep 0

compute gofr all rdf 125 1 3 1 4 2 3 3 12 3 13 12 13 3 5 5 12 #ao-ob ao-oh st-ob ob-o* ob-h* o*-h* ob-Na o*-Na 壳层数量为100，看data文件才知道1和2代表啥，这里计算的是全局尺量
fix rdf all ave/time 10 10 100 c_gofr[*] file ${filename1}.txt mode vector

 

# compute gofr all rdf 100   #计算整体rdf （global RDF），不仅仅某两个原子的原子对RDF，壳层数量为100，看data文件才知道1和2代表啥，这里计算的是全局尺量

fix rdf all ave/time 10 10 100 c_gofr[*] file ${filename1}-RDF.dat mode vector

 

 

 

run 100