# I suppose we should try to run a simulation?

units		real
atom_style	full
boundary	p p p

read_data		LAMMPS_INPUT.data &
				extra/atom/types 2 &
				extra/bond/types 1 &
				extra/angle/types 1 &
				extra/bond/per/atom 2 &
				extra/angle/per/atom 1 &
				extra/special/per/atom 2


pair_style		lj/cut/coul/cut 10.0 12.0
pair_modify		mix arithmetic
bond_style      harmonic
angle_style		harmonic
dihedral_style  fourier 

special_bonds	amber

molecule	h2o_mol H2O.txt toff 8 boff 11 aoff 21

#############################################################
include			GAFF_params.ff

pair_coeff		9 9 0.15535 3.166
pair_coeff		10 * 0.0000 0.0000

bond_coeff		12 400.0 1.0

angle_coeff		22 50.0 109.47

mass 9 15.999000
mass 10 1.008000
#############################################################

neighbor		3.0 bin
neigh_modify	delay 0 every 1 check yes

#############################################################
# Aight, lets do this already!

variable temp_seed equal 418531
variable rm_temp equal 25.0+273.15

velocity all create ${rm_temp} ${temp_seed} dist uniform mom yes rot yes units box

timestep 1.0

fix 0 all nvt temp ${rm_temp} ${rm_temp} 100.0
run 1000
unfix 0

#############################################################
reset_timestep 0

timestep 1.0

group h2o type 9 10

fix water_shake h2o shake 0.0001 50 0 b 12 a 22 mol h2o_mol
fix gcmc_nvt all nvt temp ${rm_temp} ${rm_temp} 100.0

compute_modify  thermo_temp dynamic/dof yes
compute_modify  gcmc_nvt_temp dynamic/dof yes

variable mu index -9.0
variable disp index 0.5
variable tfac equal -5.0/3.0 # (3 trans + 2 rot)/(3 trans)
fix do_gcmc h2o gcmc 10 10 0 0 418531 ${rm_temp} ${mu} ${disp} &
	mol h2o_mol pressure 0.001 tfac_insert ${tfac} group h2o shake water_shake

variable 	oxygen atom "type==9"
variable 	hydrogen atom "type==10"
group 	    oxygen dynamic all var oxygen
group 		hydrogen dynamic all var hydrogen
variable    nO equal count(oxygen)
variable    nH equal count(hydrogen)

thermo_style    custom step temp press pe ke atoms v_nO v_nH
thermo          1

run 10000