Dear All,
I have run simulation of silicon with Adri C. T. van Duin’s reaxff potential. The purpose of the simulation is to calculate thermal expansion coefficent
The problem is that the size of box is almost the same at all temperature (100k to 1000k). Silicon’s thermal expansion coefficient is positive. Therefore, the box size should increase as temperature increase.
I am not sure if it’s because this reaxff potential couldn’t predict thermal property very well or there’s something wrong with my code. There’s no research about silicon’s thermal property with this potential function. Is there anyone that has done similar simulation with this potential before ?
The potential was developed by Adri C. T. van Duin. The potential file can be found in the supporting info of paper : Integrated atomistic chemical imaging and reactive force field molecular dynamic simulations on silicon oxidation. The version of lammps I used is 11/17/2016. I have installed reaxff and kspace package.
Any reply will be appreciated.
My script is attached below.
INITIALIZATION ----------------------------
units real
dimension 3
boundary p p p
atom_style charge
atom_modify map array
ATOM DEFINITION ----------------------------
variable simu_time equal 200000
variable rest_num equal 10
variable rest_time equal {simu_time}/{rest_num}
variable itemp equal 100
variable lat equal 5.431
variable lengx equal 5*${lat}
region reg1 block 0 {lengx} 0 {lengx} 0 ${lengx} units box
create_box 1 reg1
lattice diamond ${lat}
create_atoms 1 region reg1
FORCE FIELDS ----------------------------
set type 1 charge 0
mass 1 28.086
pair_style reax/c lmp_control checkqeq yes
pair_coeff * * ffield_APL_2014_SiO Si
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
SETTINGS ----------------------------
velocity all create ${itemp} 111 dist gaussian mom yes rot yes
reset_timestep 0
timestep 0.25
variable simu_step equal round({simu_time}/dt) variable rest_step equal round({rest_time}/dt)
FIX DEFINITION ----------------------------
fix fix0 all qeq/reax 1 0.0 10.0 1.0e-6 param1.qeq
fix fix1 all npt temp {itemp} {itemp} 1000 iso 0 0 1000 tchain 10 pchain 10 nreset 100 tloop 5 ploop 5 mtk yes
variable flx equal lx
variable fly equal ly
variable flz equal lz
fix fix3 all ave/time 1 1 1 v_flx v_fly v_flzfile data.out format %20.10e
thermo 10
thermo_style custom step time pe temp lx vol pxx pyy pzz
Loop
variable i loop {rest_num}
label loop
run {rest_step}
write_restart rest_*.in
next i
jump infile_model1_equil.in loop
Regards&Thanks,
Wenting