Thank you so much for your advice. I actually want to test the stability of the BCN, so I followed your suggestion to use 3 dimension. And I also noticed that reaxff potential can only use for the system that the dimension bigger than 10 angstroms. So that, I make my system as 744350, and still give a 10k initial velocity run 10000step. The time step is 0.1 femtosecond. However, the temperature still goes to over 1000K. I visualized the trajectory, it seems like good ( see the attachment, the first step of system). Could you please help me check that where I got wrong?
Below is my new script.
Sincerely,
Peng
Heat BCN crystal, ReaxFF test
Initialization
units real
variable
variable b equal 4.26
variable a equal $b*sqrt(3.0)
variable 1_3 equal 1.0/3.0
variable 2_3 equal 2.0/3.0
variable 5_6 equal 5.0/6.0
variable 1_6 equal 1.0/6.0
Setup Box
dimension 3
boundary p p p
atom_style charge
lattice custom 1.0 &
a1 $a 0.0 0.0 &
a2 0.0 $b 0.0 &
a3 0.0 0.0 100.0 &
basis 0.0 0.0 0.0 &
basis ${1_3} 0.0 0.0 &
basis ${2_3} 0.0 0.0 &
basis 0.0 ${2_3} 0.0 &
basis {1_6} {1_6} 0.0 &
basis 0.5 ${1_6} 0.0 &
basis {5_6} {1_6} 0.0 &
basis ${1_6} 0.5 0.0 &
basis 0.5 0.5 0.0 &
basis ${5_6} 0.5 0.0 &
basis {1_3} {2_3} 0.0 &
basis {2_3} {2_3} 0.0 &
region box block 0 10 0 10 -0.25 0.25
create_box 3 box
create_atoms 1 box &
basis 2 2 &
basis 3 3 &
basis 4 2 &
basis 5 3 &
basis 7 2 &
basis 8 3 &
basis 9 2 &
basis 12 3 &
group B type 1
group N type 2
group C type 3
mass 1 10.811
mass 2 14.0067
mass 3 12.0107
Settings
pair_style reax/c NULL
pair_coeff * * ffield.reax.new B N C
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
compute reax all pair reax/c
timestep 0.1
thermo 1
thermo_style custom step temp press pe ke etotal vol density epair lx ly lz
dump 99 all custom 1 dump.reax.atom id type mass xs ys zs
Minimization
min_style cg
minimize 1.0e-10 1.0e-10 1000000 1000000
Heating Stage
velocity all create 10.0 12345 rot yes dist gaussian
fix 2 all nve
run 10000
