# Initialization
units       real
boundary    p p p
atom_style  full
neighbor    2.0 bin
neigh_modify delay 1 every 1 check yes one 20000 page 200000

# Force field styles
pair_style     lj/cut/coul/long 10.0 12.0
bond_style     harmonic
angle_style    harmonic
dihedral_style opls
kspace_style   ewald 1.0e-4
kspace_modify  gewald 0.25

# Read graphite structure
read_data data.graphite &
  extra/atom/types 3 &
  extra/bond/types 1 &
  extra/angle/types 1 &
  extra/dihedral/types 1 &
  extra/special/per/atom 6 &
  extra/bond/per/atom 4 &
  extra/angle/per/atom 6 &
  extra/dihedral/per/atom 6

# Extend box
change_box all y final -210.0 106.429025 units box

# Define regions
region slab1 block INF INF INF INF 0.0 10.0 units box
region slab2 block INF INF INF INF 50.0 60.0 units box
region wall_end block INF INF 101.0 106.429025 10.0 50.0 units box
region bulk block INF INF -100.0 -0.1 INF INF units box
region channel block INF INF INF 101 10.0 50.0 units box
region piston block INF INF -200.0 -198.0 INF INF units box

# Delete atoms in channel
delete_atoms region channel

# Create piston
lattice fcc 3.0
create_atoms 5 region piston

# Grouping
group slab1_atoms region slab1
group slab2_atoms region slab2
group wall_atoms region wall_end
group piston_atoms type 5

group fixed_slabs union slab1_atoms slab2_atoms
group fixed_wall union wall_atoms
group fixed_atoms union fixed_slabs fixed_wall

# Assign types
set group wall_atoms type 2
set group piston_atoms type 5

# Masses
mass 1 12.011
mass 2 12.011
mass 3 15.035
mass 4 14.027
mass 5 12.011

# Force field parameters
bond_coeff     1 500.0 1.54
angle_coeff    1 62.1 114.0
dihedral_coeff 1 0.0 0.0 0.0 0.0

pair_coeff 1 1 0.068443 3.407
pair_coeff 2 2 0.068443 3.407
pair_coeff 1 2 0.068443 3.407
pair_coeff 3 3 0.1947   3.75
pair_coeff 4 4 0.0910   3.95
pair_coeff 3 4 0.133    3.85
pair_coeff 1 3 0.133    3.85
pair_coeff 1 4 0.0910   3.95
pair_coeff 2 3 0.133    3.85
pair_coeff 2 4 0.0910   3.95
pair_coeff 5 5 0.068443 3.407
pair_coeff 1 5 0.068443 3.407
pair_coeff 2 5 0.068443 3.407
pair_coeff 3 5 0.133    3.85
pair_coeff 4 5 0.0910   3.95

# Freeze fixed atoms
fix freeze_slabs fixed_slabs setforce 0.0 0.0 0.0
fix freeze_wall fixed_wall setforce 0.0 0.0 0.0

# Add hexane molecules
molecule c6 C6.txt
create_atoms 0 random 100 12345 bulk mol c6 99999 overlap 0.3

# Group hexane and ensure not frozen
group hexane type 3 4
group fixed_atoms subtract fixed_atoms hexane

# Finalize groups and assign velocities
reset_atoms id
group mobile_atoms subtract all fixed_atoms piston_atoms
velocity mobile_atoms create 300.0 12345 mom yes rot yes dist gaussian loop geom
minimize 1.0e-6 1.0e-8 10000 100000

# Thermostat
compute mobTemp mobile_atoms temp
fix nvt_mobile mobile_atoms nvt temp 300.0 300.0 10.0
fix_modify nvt_mobile temp mobTemp

thermo_modify temp mobTemp
thermo 1000
thermo_style custom step temp press density lx ly lz
timestep 0.01

# Piston force
fix piston_nve piston_atoms nve
fix piston_force piston_atoms aveforce 0.0 1634.88 0.0

# Output and run
dump dmp all atom 1000 dump_graphite_channel.lammpstrj
dump_modify dmp scale no

run 200000
unfix nvt_mobile

# Save final structure
write_data data.graphite_channel_c6