# 3d LJ h2o-co2 simulation  
    #1  #c
    #2  #h
    #3  #o*
    #4  #h*
units real
dimension 3
boundary p p p
atom_style full
bond_style  harmonic 
angle_style harmonic

pair_style hybrid lj/cut/coul/long 10 9 lj/cut/tip4p/long 3 4 2 2 0.1250 10 9
kspace_style pppm/tip4p 1.0e-4

read_data  free.data

special_bonds lj/coul 0.0 0.0 0.5

pair_coeff  1    1 lj/cut/coul/long  0.06600    3.50000
pair_coeff  2    2 lj/cut/coul/long  0.03000    2.50000
pair_coeff  3    3 lj/cut/tip4p/long 0.16275    3.16435
pair_coeff  4    4 lj/cut/tip4p/long 0.00000    0.00000

pair_coeff  1    2 lj/cut/coul/long  0.04450 	3.00000 
pair_coeff  1    3 lj/cut/coul/long  0.25500 	3.03200 
pair_coeff  1    4 lj/cut/coul/long  0.00000 	1.75000 
pair_coeff  2    3 lj/cut/coul/long  0.06987 	2.83218 
pair_coeff  2    4 lj/cut/coul/long  0.00000 	1.25000 
pair_coeff  3    4 lj/cut/tip4p/long 0.00000 	0.00000 


neighbor 2 bin
neigh_modify delay 0 every 5 check yes

group C             type    1 
group H             type    2 
group OTIP          type    3 
group HYIP          type    4 

group water            type       3 4
group hydrate          molecule   3001           

compute comT all com

timestep 1.0 
thermo_style custom step temp press density pe etotal xhi yhi zhi c_comT[1] c_comT[2] c_comT[3]
thermo 1000

#energy minimization
minimize 1.0e-5 0 100000 100000

#parameter initiallized
variable ct equal 250
variable cp equal 500

#equilibrate with langevin for temperature
log   log.langevin.equilibrium

reset_timestep  0

velocity all create ${ct} 2244352 dist gaussian mom yes rot yes
velocity hydrate set 0.0 0.0 0.0

fix 0 water shake 0.00001 10 0 b 2 a 2

fix 1 all langevin ${ct} ${ct} 100.0 54324 zero yes
fix 2 all nve

fix 3 hydrate setforce 0.0 0.0 0.0

dump 1 all xyz 500 vmd1.xyz
dump_modify  1 element C H O H

run 50000

unfix 1
unfix 2
unfix 3
undump 1

#equilibrate with langevin for temperature and berendsen for pressure
log  log.npt.equilibrium
reset_timestep  0
thermo 1

velocity all zero linear

fix 4 all press/berendsen iso ${cp} ${cp} 1000.0
fix 5 all langevin ${ct} ${ct} 100.0 54324 zero yes
fix 6 all nve

dump 2 all xyz 1 vmd2.xyz
dump_modify  2 element C H O H

run 2000000

unfix 4
unfix 5
unfix 6
undump 2

#NPT with nose-hoover
log   log.npt.nose
reset_timestep  0

velocity all zero linear

fix 7 all npt temp ${ct} ${ct} 100.0 iso ${cp} ${cp} 1000.0

dump 3 all xyz 5000 vmd3.xyz
dump 4 all atom 5000 ovito3.lammpstrj
dump_modify  3 element C H O H

restart 1000000 restart.npt.*

run 100000000