variable     Ts equal 300.0
variable     Ts_end equal 700.0
variable     dt equal 2.0
variable     spring_factor equal 67.6425
variable     Tdamp equal 100*dt
variable     Pdamp equal 1000*dt
variable     n equal 10

# ---------- Initialize Simulation --------------------- 
clear 
echo screen
units        real
dimension    3
boundary     p p f
atom_style   full
#atom_modify  map array
bond_style   harmonic
angle_style  harmonic
processors   * * 1

# ---------- Create Atoms ---------------------

read_data    system.data
include      forcefield.system

# ---------- Group and Region Definition --------------------- 
group        solid type 1

group        water type 2 3
region       fixed block EDGE EDGE EDGE EDGE EDGE 1 units box
region       phan  block EDGE EDGE EDGE EDGE 1    6 units box
group        fixed region fixed
group        phan  region phan
group        wall  subtract solid fixed
group        free  subtract wall phan
group        gthermo  subtract all fixed

region       rheatflux block INF INF INF INF 6 10 units box
group        gheatflux region rheatflux

neighbor     2.0 bin 
neigh_modify delay 0 every 1 check yes exclude group fixed fixed

# ---------- Define Settings --------------------- 
compute      myPE   all pe/atom 
compute      peng water reduce sum c_myPE
compute      myKE   all ke/atom 
compute      keng water reduce sum c_myKE
variable     wateng equal "c_peng+c_keng"
compute      wtemp  water temp
compute      wtempcom  water temp/com
compute      stemp  wall  temp
compute      ptemp  phan  temp
compute      ftemp  free  temp

compute      myStress all stress/atom NULL virial
compute      flux gheatflux heat/flux myKE myPE myStress 

# ---------- Run Minimization --------------------- 

fix          spring solid spring/self ${spring_factor}
fix          freeze fixed setforce 0 0 0  #fix the atoms at two ends
fix          top_wall all wall/reflect zhi EDGE units box

thermo       1000
thermo_style custom step pe c_wtemp c_wtempcom c_stemp c_ptemp c_ftemp press pxx pyy pzz v_wateng c_flux[1] c_flux[2] c_flux[3]
min_style    cg 
minimize     1e-6 1e-8 5000 10000 

write_data   system_min.data
reset_timestep 0 

# ---------- Run Equilibration --------------------- 
velocity     gthermo create ${Ts} 13425 dist gaussian
velocity     fixed set 0.0 0.0 0.0 sum no units box
fix          fShakeSPCE water shake 0.0001 10 0 b 1 a 1

dump	     mydump1 all custom 10000 equil.*.lammpstrj id type x y z vx vy vz
dump_modify  mydump1 sort id

fix 		 loadbalance all balance 50000 1.0 shift xy 10 1.1 weight time 0.8 out balance.output

#fix          fNVEall gthermo nve
#fix          fNVTlan gthermo langevin ${Ts} ${Ts} ${Tdamp} 54687 tally yes
fix          fNVEwat water nve
fix          fNVTwat water temp/berendsen ${Ts} ${Ts} ${Tdamp}
fix          fNVEwall wall nve
fix          fNVTwall wall temp/berendsen ${Ts} ${Ts} ${Tdamp}

compute      coord  water coord/atom cutoff 6.3	2
variable     liquid_ID atom (c_coord>$n)||(c_coord==$n)
variable     vapor_ID atom (c_coord<$n)
compute      NL  water reduce sum v_liquid_ID 
compute      NV  water reduce sum v_vapor_ID

variable     vapor_ID2 atom (c_coord<12)
compute      NV2  water reduce sum v_vapor_ID2

fix		     ave_equil all ave/time 10 5000 50000 c_wtemp c_wtempcom c_stemp c_ptemp v_wateng c_peng c_NV c_NL c_NV2 c_flux[1] c_flux[2] c_flux[3] ave one file average-equil3.profile 

timestep     ${dt}
run	         500000

unfix        fNVTwat
run	         500000

unfix        fNVEwat
unfix        fNVEwall
unfix        fNVTwall
fix          fNVTpre gthermo nvt temp ${Ts} ${Ts} ${Tdamp}
run	         500000

write_data    system_equil.data
write_restart equil.restart
#reset_timestep 0

# ----------------- Main NEMD simulation: Temperature ramp in the phantom layers -----------------
undump       mydump1
dump	     mydump all custom 5000 run.*.lammpstrj id type x y z vx vy vz
dump_modify  mydump sort id

unfix        fNVTpre
#group        gnve  union free water
fix          fNVE gthermo nve
fix          tempramp phan langevin ${Ts_end} ${Ts_end} ${Tdamp} 36378 tally yes

# ----------------- Density&temperature profiles in fluid region-----------------

compute      chunk_z water chunk/atom bin/1d z lower 2.0 units box
fix          dens_z water ave/chunk 10 500 5000 chunk_z density/number density/mass c_myPE ave one file density_z.profile

compute      chunk_xz all chunk/atom bin/2d x lower 2.0 z lower 2.0 units box
fix          dens_xz water ave/chunk 10 5000 50000 chunk_xz density/mass c_myPE ave one file density_xz.profile
fix          temp_xz all ave/chunk 10 5000 50000 chunk_xz temp ave one file temp_xz.profile

fix		     ave all ave/time 10 5000 50000 c_wtemp c_wtempcom c_stemp c_ptemp v_wateng c_NV c_NL c_NV2 f_tempramp c_flux[1] c_flux[2] c_flux[3] ave one file average.profile 

compute      myRDF water rdf 100 3 3
fix          myRDF water ave/time 10000 1 10000 c_myRDF[*] file tmp.rdf mode vector

restart      1000000 run.*.restart

run          2500000

print "All done!"