#第一部分 设定模拟的基本规则
units                    real 
dimension            3
boundary             p p p
neighbor             3.0  bin
neigh_modify       every 1 delay 0 check yes
atom_style           full 
pair_style	            lj/class2/coul/long 12.5 12.5
bond_style	    class2
angle_style	    class2
dihedral_style	    class2
improper_style    class2
pair_modify	     mix sixthpower tail yes
special_bonds  	lj/coul 0 0 1.0
kspace_style           pppm 0.0001

# 第二部分  定义模拟的几何模型
read_data     jdd.data 
kspace_style           pppm 0.0001  #这是因为是三斜晶系，要再写一遍这个

# 第三部分  建立模拟物理模型Masses
mass  1   12.01115  # c1
mass  2   12.01115  # c2
mass  3   12.01115  # coh
mass  4    1.00797  # hc
mass  5    15.9994  # o
mass  6    15.9994  # oc

#第四部分 建立模拟的分组信息
# 定义上部区域（上35%）
#region    up block INF INF INF INF 76.9 115 units box
#group     up region up 
 # 定义下部区域（下35%）
#region    down block INF INF INF INF 7 45.2 units box
#group     down region down 
 # 定义中间退火区域（中间30%）
#group     anneal subtract all up down
# 定义每种原子的质量和数量变量
# 为每种原子类型创建组
group type1 type 1    # c1
group type2 type 2    # c2
group type3 type 3    # coh
group type4 type 4    # hc
group type5 type 5    # ho
group type6 type 6    # oc
group type7 type 7    # oh

# 创建包含所有原子的组
group all type 1:7

# 计算每种原子的数量和质量
variable n1 equal count(type1)
variable n2 equal count(type2)
variable n3 equal count(type3)
variable n4 equal count(type4)
variable n5 equal count(type5)
variable n6 equal count(type6)
variable n7 equal count(type7)
# 定义阿伏伽德罗常数
variable NA equal 6.02214076e23

# 计算总质量 (g/mol)
variable total_mass equal 12.01115*v_n1 + 12.01115*v_n2 + 12.01115*v_n3 + 1.00797*v_n4 + 1.00797*v_n5 + 15.9994*v_n6 + 15.9994*v_n7

# 体积转换
variable vol_A3 equal vol                           # 体积（ų）
variable vol_cm3 equal v_vol_A3*1.0e-24            # 转换为 cm³

# 计算正确的密度（考虑阿伏伽德罗常数）
variable calculated_density equal ${total_mass}/(v_NA*v_vol_cm3)

# 输出信息
print "Volume = ${vol_A3} A^3"
print "Volume = ${vol_cm3} cm^3"
print "Total mass = ${total_mass} g/mol"
print "Calculated density = ${calculated_density} g/cm^3"
print "LAMMPS density = $(density) g/cm^3"

quit

#第五部分 弛豫
thermo 1000
thermo_style custom step temp press vol density pe ke etotal 
thermo_modify flush yes
dump 1 all custom 10000 dump.lammpstrj id mol type mass x y z

#能量最小化
minimize 0.0 1.0e-8 1000 100000

#初始温度
compute myT all temp 
velocity all create 300 262027 rot yes dist gaussian

# 设置时间步长
timestep 1

#设置NVT弛豫

fix 1 all nvt temp 300 300 100 

#开始循环模拟
label loop
variable i loop 1000
run 1000
write_data temp.data
next i
jump SELF loop


#数据输出
write_data  chiyuu_structure.data
write_restart restartt.equil