######################################################################################
variable	t equal 353
variable 	p equal 80
variable	phi equal 0.875
variable	pp equal 2000
variable	chemu equal -2.81154

variable        height1 equal 7.7001
variable        height2 equal 51.6999

######################################################################################

variable 	seed equal 31798.0000 

units		metal
atom_style      full
boundary 	p p f
atom_modify	map array

pair_style	lj/cut 14.0
pair_modify 	mix arithmetic

######################################################################################
read_data	4.4nm_full.txt
molecule	CH4 CH4.txt

group 		wall type 1
group		gas type 2
######################################################################################
variable        rxlo equal xlo+0.0001
variable        rxhi equal xhi-0.0001
variable        rylo equal ylo+0.0001
variable        ryhi equal yhi-0.0001
variable        rzlo equal zlo-0.0001
variable        rzhi equal zhi+0.0001

region		1 block ${rxlo} ${rxhi} ${rylo} ${ryhi} ${rzlo} ${height1} units box side in
group		wallbottom region 1
region          2 block ${rxlo} ${rxhi} ${rylo} ${ryhi} ${height2} ${rzhi} units box side in 
group		wallup region 2
region          active block ${rxlo} ${rxhi} ${rylo} ${ryhi} ${height1} ${height2} units box side in

group		wall union wallbottom wallup
set		group wall type 1
group		gas subtract all wall
set		group gas type 2

neighbor	1.5 bin 
neigh_modify	every 5 delay 0

compute         allperatom gas stress/atom NULL
compute         allp gas reduce sum c_allperatom[1] c_allperatom[2] c_allperatom[3]
variable        allpress equal -(c_allp[1]+c_allp[2]+c_allp[3])/(3*vol)

compute  mdtemp gas temp
##### Initial Configuration Stabilization ########################################################
velocity        gas create $t ${seed} mom yes rot yes dist gaussian

fix             setforce wall setforce 0.0 0.0 0.0
minimize        1.0e-3 1.0e-3 10000 10000

neighbor	1.5 bin
neigh_modify	every 1 delay 0
reset_timestep	0
##### GCMC #################################################
#fix            ID group-ID gcmc N X M type seed T mu displace keyword values 
fix             GC gas gcmc 100 300 300 0 ${seed} $t ${chemu} 0.1 mol CH4 pressure $p fugacity_coeff ${phi} group gas region active full_energy tfac_insert 1.00 maxangle 180
compute_modify	mdtemp dynamic yes

variable 	nm equal count(gas)
variable 	totene equal etotal 

fix             check gas ave/time 1 10 300 c_mdtemp v_allpress v_totene v_nm ave one file gcmccheck.txt
fix             check2 gas ave/time 1 5000 10000 c_mdtemp v_allpress v_totene v_nm ave one file ave.txt

thermo		1000
thermo_style	custom step v_nm c_mdtemp v_allpress etotal
thermo_modify	format float %12.10g

timestep	0.001
run		100000