#Using fix gcmc combined with fix nvt to add spc water in empty box
units		real
dimension 3
atom_style	full
boundary    p p p
region myregion block 0 80  0  80  0  80 units box
create_box 3 myregion bond/types 1 angle/types 1 extra/bond/per/atom 2 &
             extra/angle/per/atom 1 extra/special/per/atom 2
pair_style  lj/cut/coul/cut  20.0
#Using long range coulomb causes problems!
#pair_style  lj/cut/coul/long/omp  20.0
#kspace_style ewald/omp 1.0e-4
pair_modify mix arithmetic
bond_style   harmonic
angle_style  harmonic

molecule h2o mol.water.txt #toff 5 #spc water
mass    1     15.999   # OW
mass    2	     1.008   # HW1
mass    3      1.008   # HW2
pair_coeff    1    1     0.1554000       3.16557       20.0    # OW
pair_coeff    2    2     0.00000e+00     0.00000e+00   20.0    # HW1
pair_coeff    3    3     0.00000e+00     0.00000e+00   20.0    # HW2
bond_coeff  	1   450    1.000  #spc water
angle_coeff   1    55  109.470  #spc water

group water type 1:3
fix  wshake  water shake 0.0001 50 0 b 1 a 1 mol h2o
variable temp equal 375.0
#? shake keyword in fix gcmc causes LAMMPS getting stuck in computing !
#  remove shake keyword in fix gcmc works fine !
fix wgcmc water gcmc  1000  1000 1000 0 2949499 ${temp} -0.48 0.5 &
    mol h2o maxangle 180 #shake wshake 
fix wnvt  water nvt temp ${temp} ${temp} 100

compute T4 water temp/com
compute_modify T4 dynamic/dof yes
compute_modify thermo_temp dynamic/dof yes
thermo 100
thermo_style  custom step atoms temp press density evdwl ecoul c_T4

dump 5 all xyz  1000 ./water.*.xyz
dump_modify  5 element OW HW1 HW2
dump_modify  5 format line "%5s %8.3f  %8.3f  %8.3f"

timestep 1.0
run 500000