#-----------------Basic setup-------------
clear
units           real
dimension       3
boundary        s s s    
                   
 atom_style    full
 #bond_style      class2
 #bond_coeff     1 1.0 100.0 100.0 80.0
 #dihedral_style  harmonic
 #improper_style  harmonic

#-------------Potential-------------------

pair_style reax/c NULL 
read_data      m2.data
#pair_coeff * * ffield.reax C C O N H H





#---------------Structure-----------------

region         top block INF INF INF INF INF -50 units box  
region         bottom block INF INF INF INF 50 INF units box
group          g_cnt type 1
group          g_top region top
group          g_bottom region bottom
group          g_boundary union g_top g_bottom 
group          g_body subtract all g_boundary



#------------NPT relaxation at 300K-----------       
reset_timestep 0
timestep 0.001
velocity         g_body create 300 4952925 dist gaussian rot yes units box
fix              1 all npt temp 300 300 0.01 z 1 1 0.1
thermo 100
thermo_style    custom step temp c_body_temp pe lx ly lz etotal pzz c_strAll vol 
thermo_modify   lost ignore norm yes

variable length equal "lz"
print "Lattice constant (Angstoms) = ${length};"


fix              3 g_boundary setforce 0.0 0.0 NULL
dump             4 all custom 500 trajectory_s1.lammpstrj id x y z ix iy iz
dump_modify      4 element C

run              10000