# 用速度重标法计算石墨烯热导率的LAMMPS脚本

# -------------定义系统单位-------------
units       metal       # 采用metal单位系统（记住：能量eV；长度：A；时间：ps）
boundary    p p p       # 采用周期性边界条件
atom_style  atomic      # 采用原子类型
neighbor    0.3 bin     # 邻居半径
neigh_modify every 1 delay 10 check yes  #设置邻居列表的更新方式

# 定义几个参数
variable  T equal 300 # 温度为 60 K
variable  A equal 1.42 # 晶格常数大概为 1.42 A
variable  DT equal 0.0002 # 积分步长为 0.2 fs
variable  POWER equal 0.005 # 与热流对应的功率 0.05 eV/ps
# -------------建模--------------------
# a1 a2 a3代表x y z三个方向上晶格的间距（单位长度）
lattice custom  ${A} a1  3  0  0   a2   0 1.732  0 a3 0 0 2.357   &
basis   0     0   0  & 
basis   0.333 0   0  & 
basis   0.5   0.5 0  &
basis   0.833 0.5 0

region graphenesheet block 0 255.78 0 24.12 -0.5 0.5 units box   #定义一个名为graphenesheet的区域，大小为50*50*1个晶格,一共50*50*1*2=5000个原子
create_box      1 graphenesheet                           #在graphenesheet区域内先造一个盒子，原子类型为1
create_atoms    1 region graphenesheet                    #在盒子内填满原子
mass            1 12.0107                                      #设置原子1的质量为12
change_box      all z delta -150 150 units box              #填充真空层

# -------------设置势函数--------------
pair_style  tersoff
pair_coeff  * * SiC.tersoff C

write_data      1Model.data
#-------------速度初始化、能量最小化---------------
timestep    ${DT}  # 设置积分步长
velocity    all create ${T} 89684 dist gaussian #设置原子的初始速度，温度为300K，随机种子为12345
min_style   cg
minimize    1e-10 1e-10 1000 1000  

write_data      2min.data

#-------------在NVT下进行弛豫-----------
dump        1 all custom 10000 4graphene.lammpsstrj id type x y z fx fy fz # 输出轨迹文件graphene.lammpsstrj
fix         NVT all nvt temp ${T} ${T} 0.02 # Nose-Hoover 热浴
thermo_style custom step temp press pe ke lx ly lz density # 观察温度和压强的变化
thermo      100000 # 1000 步输出一次

# 根据以上设置跑若干步，目的是使系统达到热力学平衡态
run 1000000  # 跑 20000 步
undump 1
write_data  3chiyu.data

unfix NVT   # 撤掉之前的 NVT 系综
fix  NVE all nve  # 换上 NVE 系综
run 100000

reset_timestep 0

# 定义热源（heat source）和热汇（heat sink）
region SOURCE block 0 1 INF INF INF INF units box # 第1个block是热源
region SINK block 25 26 INF INF INF INF units box # 第11个block是热汇

# 在热源注入热量，在热汇导出热量，使系统达到非平衡稳态

fix HEAT_IN  all ehex 1 ${POWER} region SOURCE # 注入热量
fix HEAT_OUT all ehex 1 -${POWER} region SINK # 导出热量
# 可将ehex换成heat。使用前者能量守恒满足得更好。

# 计算单原子“温度”
compute  KE all ke/atom    # 单原子动能
variable KB equal 8.617e-5 # 玻尔兹曼常数，eV/ps
variable TEMP atom c_KE/1.5/${KB} # 单原子温度

# 下面两个命令是联合起来计算块温度（block temperatures）的。
# 首先定义一个计算，将体系在x方向平均分为1/0.02=50块，
# 然后计算每个块的温度的时间平均，将结果记录到文件。
compute BLOCKS all chunk/atom bin/1d x lower 0.02 units reduced
fix     T_PROFILE all ave/chunk 10 10000 100000 BLOCKS v_TEMP file temp${T}.txt
thermo      100000        # 每100000步输出一次热力学信息
thermo_style custom step temp press lx ly vol
# 输出自定义信息，步数、温度、压强、三个方向的热导率、三个方向的热导率的平均值
dump        1 all custom 100000 5conductivity.lammpsstrj id type x y z fx fy fz # 输出轨迹文件graphene.lammpsstrj
# 在施加热流的情况下跑若干步
run 5000000             #输出20次