#======General Set=======
units		    metal                                 #Note the Units!!
boundary	    p p p
atom_style	   full                    #Full to include bonds and angles
molecule   water   wateroffset.molecule # molecule templates
#===========Variable defination===============
variable        slit_thickness   equal   10.0  # Must be larger than 3.0A
variable        lattice_z        equal   5.4052
variable        slit_zlo         equal  2.0*${lattice_z}
variable        slit_zhi         equal  ${slit_zlo}+${slit_thickness}
variable        box_zhi         equal  2*${slit_zlo}+${slit_thickness}
variable        box_xhi         equal  4*7.3695
variable        box_yhi         equal  4*4.25478
variable        spring_up         equal  ${slit_zhi}+0.5


variable        chem         equal  -35.0/96.15 #-65.5/96.15
#======Crystal structure of SiO2_quartz=======
lattice 	custom  4.913   origin 0.0 0.0 0.0 &
		a1 1 0 0  a2 -0.5 0.8660254 0  a3 0 0 1.100183  &      
		basis      0.53   0.53   0.000  &
		basis      0.47   0.00   0.333  &
		basis      0.0    0.47   0.6667 &
 		basis      0.414  0.268  0.214  &
		basis      0.586  0.854  0.119  &
		basis      0.146  0.732  0.453  &
		basis      0.732  0.146  0.547  &
		basis      0.268  0.414  0.786  &
		basis      0.854  0.586  0.881  
    
#======Region Defination=======    
region          total    prism           0 ${box_xhi} 0 ${box_yhi} 0 ${box_zhi} 0 0 0  units box 
#======Create simulation box=======    
create_box   5 total bond/types 1 angle/types 1  extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2

region          up       prism      0 ${box_xhi} 0 ${box_yhi} ${slit_zhi} ${box_zhi}   0 0 0 units box# upper layer
region          down       prism      0 ${box_xhi} 0 ${box_yhi} 0 ${slit_zlo}   0 0 0 units box# down layer 


#======Create silica layers=======    
create_atoms   2          region   up  basis  1 1  basis 2 1 basis 3 1 &
basis  4 2  basis 5 2 basis 6 2 &
basis  7 2  basis 8 2 basis 9 2
create_atoms   2          region   down  basis  1 1  basis 2 1 basis 3 1 &
basis  4 2  basis 5 2 basis 6 2 &
basis  7 2  basis 8 2 basis 9 2 

mass          1          28
mass          2          16
mass          3          1
mass          4          16
mass          5          1

group water_group type 4 5
delete_atoms group water_group
variable        chem         equal  -35.0/96.15 #-65.5/96.15
variable        N_water       equal     count(water_group)/3







#======Change box type======= 
#change_box all ortho remap units box











#======bond and angle defination for water======= 
kspace_style    ewald 1.0e-5
#======pair defination for SiO2 and LJ between water molecules======
pair_style hybrid tersoff lj/cut 10.0 lj/cut/coul/long   10.0
pair_coeff     * * tersoff SiO.tersoff Si O  NULL NULL NULL 
pair_coeff      1         3       lj/cut     0.005107539    2.1132    
pair_coeff      2         3       lj/cut       0.003636769         1.77765       
pair_coeff      3         3       lj/cut     0.002         0.4
pair_coeff      1         4       lj/cut/coul/long      2.3E-5     3.23375
pair_coeff      1         5       lj/cut/coul/long      0.0    0.0   
pair_coeff      2         4       lj/cut/coul/long      0.0068      3.1655  
pair_coeff      2         5       lj/cut/coul/long       0.0         0.0    
pair_coeff      4         4       lj/cut/coul/long     0.0068       3.1655  
pair_coeff      4         5       lj/cut/coul/long    0.0         0.0
pair_coeff      5         5       lj/cut/coul/long     0.0         0.0
pair_coeff      3         4       lj/cut/coul/long      0.0    0.0 
pair_coeff      3         5       lj/cut/coul/long     0.0    0.0 






#=====Neighbor list===============
neighbor       2.0 bin
neigh_modify    check yes

timestep	        0.001
dump pre all custom 100 ./dump/sheartest.* id type x y z 


#=====================Fix GCMC====================================================================
#==========fix gcmc==========
fix       3  water_group gcmc 1 1 1 0 34565 300.0 ${chem} 1.0 mol water maxangle 18.0 
thermo 1
thermo_style custom  step v_N_water lx ly lz temp pe etotal pzz f_3[1] f_3[2] f_3[3] f_3[4] f_3[5] f_3[6] f_3[7] f_3[8]
run 10000 #Need to adjust