units             real
atom_style        full
boundary          f p p

read_data         water.data

pair_style        reax/c NULL checkqeq yes
pair_coeff        * * ffield H O 

timestep          0.5
thermo_style      custom etotal pe ke temp evdwl press pxx pyy pzz pxy pxz pyz lx ly lz
thermo_modify     line multi flush yes
thermo            10

variable          ti equal 300.0

dump              1 all custom 500 traj.dump id type x y z q

fix               0 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
reset_timestep    0

# Barostat
variable          xlo equal -4.0 
variable          xhi equal 24.0 
print             "xdim ${xlo} ${xhi}"
variable          varly equal ly
variable          varlz equal lz

# External Pressures
variable          Prext equal 1000.0
variable          Plext equal 1000.0

# Wall dynamics
variable          v0l equal 0.005
variable          v0r equal -1.*${v0l} 
variable          dv equal 0.001
variable          vl equal ${v0l}
variable          wl equal ${xlo}
variable          vr equal ${v0r} 
variable          wr equal ${xhi}

# Wall
variable          epsilon equal 1.0 
variable          sigma equal 1.0
variable          cutoff equal 2.5

variable          timer loop 1 10000000

variable          nrun equal 10
variable          halfnrun equal ${nrun}/2.

fix               1 all npt temp ${ti} ${ti} 100 y 1.0 1.0 1000 z 1.0 1.0 1000 
fix               2 all wall/harmonic xlo ${wl} ${epsilon} ${sigma} ${cutoff} xhi ${wr} ${epsilon} ${sigma} ${cutoff} units box 
fix               4 all ave/time 1 ${halfnrun} ${nrun} f_2[1] f_2[2] v_varly v_varlz file fprofile.profile
restart           50000 restart.*.rst

# Equilibration
run               ${nrun} 

variable          NA2toAtm equal 68568.415 
variable          aa equal f_4[3]*f_4[4]
variable          Pleft equal ${NA2toAtm}*abs(f_4[1])/${aa}
variable          Pright equal ${NA2toAtm}*abs(f_4[2])/${aa}
variable          gl equal (${Plext}-(${NA2toAtm}*abs(f_4[1])/${aa}))/${Plext}
variable          gr equal (${Prext}-(${NA2toAtm}*abs(f_4[2])/${aa}))/${Prext}

fix               7 all print ${nrun} ">>> Pleft ${Pleft} / ${Plext} Pright ${Pright} / ${Prext} Area ${aa}"
fix               8 all print ${nrun} ">>> Gleft ${gl} Vleft ${vl} Wleft ${wl} Gright ${gr} Vright ${vr} Wright ${wr}"

label             loop

variable          vl delete
variable          vl equal ${v0l}+${dv}*${gl}
variable          vr delete
variable          vr equal ${v0r}-${dv}*${gr}

variable          wl delete
variable          wr delete
variable          wl equal ${xlo}+${vl}*${timer}
variable          wr equal ${xhi}+${vr}*${timer}
next              timer
variable          v0l delete
variable          v0r delete
variable          v0l equal ${vl}
variable          v0r equal ${vr}

unfix             2
fix               2 all wall/harmonic xlo ${wl} ${epsilon} ${sigma} ${cutoff} xhi ${wr} ${epsilon} ${sigma} ${cutoff} units box 
run               ${nrun}   

jump              relax.in loop